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DESCRIPTION 

ULTRASON IC DIAGNOSIS SYSTEM AMH ctTPAT N DISTRIBUTION DISP LAX 

METHOD V 



5 Technical Field 

The present invention relates to an ultrasonic 
diagnosis system and a strain distribution display method 
which allow the user to make quantitative measurement of the 
rigidity of the tissue using an ultrasonic diagnosis 
10 apparatus. 

Background Art 

Ultrasonic diagnosis is applied not only to observation 
of the tissue structure but also to the field (ultrasonic 
tissue characterization) where physical quantities within 

15 the tissue such as the sound speed, the damping coefficient, 
and so forth, are measured and furthermore, diagnosis images 
are created based upon the physical quantities thus measured. 
As a part of the field, the field is known wherein the 
rigidity of the tissue, i.e., the elastic property is 

20 measured. The aforementioned field is being intensely 

studied since the elastic property of the tissue has close 
relation to the pathological situation. For example, it is 
known that the tissue affected by: sclerosing tumors such as 
mammary cancer, thyroid cancer, and so forth; liver 

25 cirrhosis; arterial sclerosis; and so forth, exhibits 
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greater rigidity than the normal tissue. Conventionally, 
the rigidity of the tissue is detected by touch. However, 
detection by touch has the disadvantage of difficulty in 
objective analysis, requires skill of the surgeon, and has 
5 the limitations that only the affected tissue having a 

certain size or more and positioned near the body surface 
can be detected. 

On the other hand, a method is known wherein static 
pressure is applied to the body surface so as to compress 

10 and deform the tissue, and the strain within the tissue 
corresponding to the applied pressure is measured using 
ultrasonic wave in order to estimate the elastic property of 
the tissue (J. Ophir, I. Cespedes, H. Ponnekanati, Y. Yazdi, 
and X. Li, "Elastography : A quantitative method for imaging 

15 the elasticity of biological tissue". Ultrasonic Imaging, 

Vol. 13, pp. 111-134, 1991). The conventional technique has 
been developed based upon the fact that the tissue having 
great rigidity exhibits small strain thereof under pressure, 
and on the other hand, the tissue having small rigidity 

20 exhibits great strain thereof under pressure. That is to 
say, with the aforementioned conventional method, static 
pressure is applied to the tissue, and the elastic property 
of the tissue is estimated based upon the strain 
distribution within the tissue under pressure thus applied. 

25 Specifically, normal measurement of ultrasonic echo 
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signals (RF signals without application of pressure) is made 
using an ultrasonic diagnosis apparatus with an ultrasonic 
probe without pressure applied to the tissue through the 
ultrasonic probe. Subsequently, the surgeon applies 
5 pressure to the tissue through the ultrasonic probe to a 

slight degree (around several percent), following which the 
ultrasonic echo signals (RF signals under pressure) passing 
through the tissue to which pressure is applied are measured. 
Then, the displacement distribution which represents 

10 displacement of each point of the tissue due to the pressure 
thus applied is estimated based upon the RF signals with and 
without application of pressure of the tissue using the 
spatial correlation method. 

The spatial correlation method has a mechanism wherein 

15 the displacement distribution within the tissue under the 

applied pressure is estimated based upon the RF signals (or 
envelope signals of the RF signals) with and without 
application of pressure applied to the tissue by template 
matching using a two-dimensional correlation function. That 

20 is to say, a two-dimensional correlation window (template) 
having a certain size is applied to the RF signal data 
corresponding to the tomographic data without pressure 
applied to the tissue so as to estimate displacement of a 
desired measurement point on the two-dimensional surface by 

25 detecting the maximum correlation between the RF signal data 
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to which the correlation window has been applied and the RF 
signal data under pressure applied to the tissue using the 
autocorrelation processing. The aforementioned 
autocorrelation processing is performed for each measurement 
5 point set in the shape of a grid, for example, whereby the 
strain distribution is estimated. In general, the 
processing using the spatial correlation method is performed 
with poorer precision of displacement detection in the 
horizontal direction (scanning direction of the ultrasonic 

10 beam) than in the axial direction due to rougher sampling in 
the horizontal direction than in the axial direction. As 
described above, the spatial correlation method has the 
advantage of enabling the user to estimate the two- 
dimensional displacement vector. Furthermore, while the 

15 aforementioned spatial correlation method has the 

disadvantage of precision of the estimated displacement 
being limited by the sampling pitch, the spatial correlation 
method has the advantage of enabling the user to estimate 
the displacement distribution even in a case wherein the 

20 tissue is greatly deformed (e.g., around 5%). However, the 
spatial correlation method has the disadvantage of being 
calculation-intensive for the spatial correlation processing, 
leading to difficulty in processing in real time, unlike the 
conventional ultrasonic diagnosis. 

25 Accordingly, it is an object of the present invention 



- 5 - 



to provide a method for obtaining the displacement 
distribution, strain distribution, and elastic modulus 
distribution, in real time. 

Disclosure of Invention 
5 In order to solve the aforementioned problems, an 

ultrasonic diagnosis system according to the present 
invention for obtaining the displacement of the tissue of 
the subject based upon the reflected echo signals (RF 
signals) by measurement of the subject with and without 

10 application of pressure using an ultrasonic probe comprises: 
storage means for storing the properties of signals such as 
the envelope signals thereof detected with said ultrasonic 
probe; correlation computing means for calculating the 
correlation coefficient between said properties with and 

15 without pressure applied to the subject and the phase 

difference between the received signals with and without 
application of pressure, based upon said properties stored 
in said storage means with and without pressure applied to 
said subject; and displacement computing means for 

20 calculating the displacement of each measurement point due 
to said application of pressure based upon the correlation 
coefficient and the phase difference between the RF signals 
with and without application of pressure thus obtained by 
the correlation computing means. Furthermore, the 

25 ultrasonic diagnosis system according to the present 
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invention may include strain computing means for calculating 
the strain distribution of the tissue of the subject by 
making spatial differentiation of the displacement at each 
measurement point, and display means for displaying the 
5 strain distribution thus obtained. 

As descried above, with the ultrasonic diagnosis system 
according to the present invention, the displacement of each 
measurement point is calculated based upon the correlation 
between the properties such as the envelope signals with and 

10 without application of pressure, thereby enabling real-time 
estimation of the displacement distribution. Furthermore, 
with the ultrasonic diagnosis system according to the 
present invention, the position of each measurement point 
may be obtained so as to exhibit the maximum correlation 

15 coefficient between the envelope signals with and without 
application of pressure for the measurement points by 
varying said measurement points in said ultrasonic beam 
direction at a pitch half the wavelength of said ultrasonic 
signals, thereby solving a problem of aliasing which is the 

20 disadvantage in the Doppler method. 

Note that the correlation computation wherein the 
position of each measurement point is calculated so as to 
exhibit the maximum correlation coefficient between the 
envelope signals with and without application of pressure 

25 for the measurement points by varying the phase of the 
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envelope signals with application of pressure along the time 
axis so as to obtain the autocorrelation function of the 
envelope signals with application of pressure for each 
measurement point leads to great calculation time, often 
5 leading to a problem that real-time calculation cannot be 
made. 

Accordingly, with the ultrasonic diagnosis system 
according to the present invention, first, the 
autocorrelation functions of the envelope signals with and 

10 without application of pressure are preferably calculated. 

Then, the correlation coefficient is obtained by varying the 
phase difference between the autocorrelation functions thus 
obtained corresponding to the predetermined variation of the 
measurement points, e.g., by varying the phase difference 

15 thereof at a pitch half the wavelength of the ultrasonic 
signals. This reduces calculation time for displacement 
computation, thereby enabling high-speed processing 
Furthermore, the ultrasonic diagnosis system according to 
the present invention may further include: storage means for 

20 storing the envelope signals of the quadrature -detected RF 
signals; correlation computing means for calculating the 
position of each measurement point which exhibits the 
maximum correlation coefficient between the envelope signals 
with and without application of pressure for the measurement 

25 points surrounded by a two-dimensional correlation window; 
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and displacement computing means for obtaining at least two- 
dimensional displacement of each measurement point due to 
application of pressure based upon the position of each 
measurement point which exhibits the correlation coefficient 
5 and the phase difference thus obtained by the correlation 
computing means . 

That is to say, the method which is referred to as 
"Combined Autocorrelation method (CA method)" is proposed in 
the present specification. As described above, the Combined 

10 Autocorrelation method has the advantages of two-dimensional 
and three-dimensional displacement measurement in the 
spatial correlation method with a correlation window, and 
the advantage of real-time and high-precision calculation in 
the Doppler method. The Combined Autocorrelation method 

15 allows the user to estimate the displacement distribution 
regardless of the displacement in the horizontal direction 
to a certain degree. In this case, the two-dimensional 
directions may comprise the ultrasonic-beam direction where 
the ultrasonic signals are received with the ultrasonic 

20 probe and the ultrasonic-beeim scanning direction. 

Furthermore, with the ultrasonic diagnosis system according 
to the present invention, the position of each measurement 
point which exhibits the maximum correlation coefficient is 
preferably obtained by varying the measurement points in the 

25 ultrasonic-beam direction at a pitch of half the wavelength 
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of the ultrasonic signals, and in the ultrasonic-beam 
scanning direction at the ultrasonic-beam pitch. Note that 
while the pitch at which the measurement points are varied 
in the ultrasonic-beam direction according to the present 
5 invention is not restricted to half of the wavelength of the 
ultrasonic signals, a pitch smaller than half of the 
wavelength of the ultrasonic signals is preferably employed. 

In order to further improve calculation speed of the 
displacement computation, first, the autocorrelation 

10 functions of the envelope signals with and without 

application of pressure are preferably calculated. Then, 
the position of each measurement point which exhibits the 
maximum correlation coefficient is obtained by varying the 
phase difference between the autocorrelation functions in a 

15 range corresponding to the predetermined variation of the 
measurement points . 

Furthermore, the displacement computation according to 
the present invention may be applied not only to two- 
dimensional calculation, but also to three-dimensional 

20 calculation. With an three-dimensional arrangement 

employing an ultrasonic probe having a one -dimensional array 
structure, the frame data stored in the storage means 
comprises volume data formed of multiple frame data sets 
each of which serves as slice data. On the other hand, with 

25 an three-dimensional arrangement employing an ultrasonic 
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probe having a two-dimensional array structure, the data 
contains the envelope signals obtained by scanning in the 
slice direction. The correlation computing means obtain the 
position of each measurement point which exhibits the 
5 maximum correlation coefficient between the envelope signals 
with and without application of pressure for the measurement 
points surrounded by a three-dimensional correlation window 
by varying the measurement points surrounded by the three- 
dimensional correlation in three-dimensional directions as 

10 to the volume data. At the same time, the correlation 

computing means calculate the phase difference between the 
RF signals with and without application of pressure. In 
this case, the three-dimensional directions may comprise the 
ultrasonic-beam direction where the ultrasonic signals are 

15 received with the ultrasonic probe, the ultrasonic-beam 

scanning direction, and the slice direction orthogonal to 
the aforementioned two directions. Furthermore, the 
correlation computing means preferably obtain the phase 
difference between the RF signals with and without 

20 application of pressure in the ultrasonic-beam direction, in 
the ultrasonic-beam scanning direction, and in the slice 
direction orthogonal to the aforementioned two directions . 
Furthermore, the high-speed processing method described 
above may be applied to an three-dimensional arrangement. 

25 Furthermore, the calculation described above may be made by 
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varying the measurement points in the slice direction in 
units of the slice pitch of the ultrasonic beam. 

Furthermore, an method for obtaining the elastic 
modulus distribution according to the present invention may 
5 include elastic modulus computation means for creating at 
least a two-dimensional or three-dimensional finite element 
model by dividing the subject into a finite number of 
elements, and computing the elastic modulus distribution 
based upon the information used for creating the model and 

10 the strain distribution thus obtained. Furthermore, the 

elastic modulus distribution thus obtained may be displayed 
with the display means. In this case, the elastic modulus 
computing means preferably create a three-dimensional finite 
element model by dividing the tissue of the subject into a 

15 finite number of rectangular parallelepiped elements on the 
assumption that the tissue of the subject exhibits isotropic 
elasticity and near-incompressibility . Furthermore, the 
elastic modulus computing means preferably compute the 
elastic modulus distribution based upon the information 

20 regarding the aforementioned strain distribution using the 
elastic equation on the assumption that the tissue of the 
subject exhibits the uniform elastic modulus, the uniform 
stress, and the uniform strain, for each element. 

As described above, calculation according to the 

25 present invention is made on the assumption that the tissue 
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exhibits isotropic elasticity. The reason is that the 
relation between the stress and strain exhibits near- 
linearity under static external pressure applied to the 
tissue. Accordingly, approximate calculation can be made on 
5 the assumption that the tissue serves as an elastic model. 
In addition, the tissue exhibits isotropic properties, and 
accordingly, calculation according to the present invention 
may be made on the assumption that the tissue serves as an 
isotropic elastic model. On the other hand, with the 

10 present invention, calculation is made on the assumption 

that the tissue exhibits near-incompressibility . The reason 
is that if calculation is made on the assumption that the 
tissue exhibits the complete incompressibility , i.e., 
calculation is made with the constant Poisson's ratio of 0.5 

15 within the tissue, the elastic equation becomes a special 
case, leading to a problem that calculation cannot be made 
using the finite element method. Note that with the present 
invention, calculation may be made with a uniform Poisson's 
ratio within the tissue. In this case, only the Young's 

20 modulus should be estimated for estimation of the elastic 

modulus distribution, thereby reducing the inverse problem. 
Note that the Poisson's ratio exhibits sufficient uniformity 
within the tissue, and accordingly, with the present 
invention, calculation is preferably made with the Poisson's 

25 ratio of 0.49. The elastic modulus distribution computation 
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according to the present invention enables reconstruction of 
the elastic modulus distribution based upon the strain 
distribution in the axial direction alone which can be 
computed with high precision, thereby enabling stable 
5 computation of the elastic modulus distribution. 

Brief Description of the Drawings 
Fig. 1 is a diagram for describing the mechanism of an 
ultrasonic diagnosis apparatus . 
10 Fig. 2 is a diagram which shows a specific example of a 

tissue elasticity measurement method by application of 
static pressure, and the mechanism of the tissue elasticity 
measurement method by application of static pressure. 

Fig. 3 is a diagram which shows the mechanism of the 
15 spatial correlation method. 

Fig. 4 is a diagram which shows the mechanism of the 
Doppler method. 

Fig. 5 is a diagram which shows the mechanism of the 
combined autocorrelation method. 
20 Fig. 6 is a block diagram which shows a circuit 

configuration for executing basic algorithm of the combined 
autocorrelation method. 

Fig. 7 is a block diagram which shows the schematic 
configuration of an ultrasonic diagnosis system according to 
25 an embodiment of the present invention. 
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Fig. 8 is a flowchart which shows basic algorithm of 
the three-dimensional combined autocorrelation method. 

Fig. 9 is a flowchart which shows the basic algorithm 
of the three-dimensional combined autocorrelation method 
5 employed in the ultrasonic diagnosis system according to the 
present invention, and a part of the processing shown in Fig. 
7 in detail. 

Fig. 10 is a flowchart for making detailed description 
regarding the combined autocorrelation method with improved 
10 calculation speed, which is performed in Step S15 shown in 
Fig. 9. 

Fig. 11 is a block diagram which shows a circuit 
configuration for executing the basic algorithm of the 
three-dimensional combined autocorrelation method employed 
15 in the ultrasonic diagnosis system according to the present 
invention . 

Fig. 12 is a schematic diagreim which shows the 
procedure of simulation. 

Fig. 13 is a diagram which shows an example of a point 
20 spread function for each point with the ultrasonic center 
frequency of 5 . 0 MHz . 

Fig. 14 is a schematic diagram which shows a tissue 
model . 

Fig. 15 is a diagram which shows the error of estimated 
25 strain with each estimating method due to displacement in 
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the horizontal direction. 

Fig. 16 is a diagram which shows the strain 
distribution estimated with each estimating method (combined 
autocorrelation method, expanded combined autocorrelation 
5 method, spatial correlation method) in a case of 

displacement of 0.0 mm in the horizontal direction. 

Fig. 17 is a diagram which shows the strain 
distribution estimated with each estimating method (combined 
autocorrelation method, expanded combined autocorrelation 
10 method, spatial correlation method) in a case of 

displacement of 0.4 mm in the horizontal direction. 

Fig. 18 is a schematic diagram which shows a tissue 
model used for simulation of compression in a slant 
direction. 

15 Fig. 19 shows estimated results of the strain 

distribution obtained from simulation of simple compression 
of the tissue model in the axial direction. 

Fig. 20 shows estimated results of the strain 
distribution obtained from simulation of compression of the 
20 tissue model in a slant direction. 

Fig. 21 is a diagram which shows two tissue model 
examples each of which have a three-dimensional structure. 

Fig. 22 is a first diagram which shows estimated 
results with an inclusion-containing model. 
25 Fig. 23 is a second diagram which shows estimated 
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results with the inclusion-containing model. 

Fig. 24 is a first diagram which shows estimated 
results with a layer- structure model. 

Fig. 25 is a second diagram which shows estimated 
5 results with the layer- structure model. 

Fig. 26 is a diagram which shows a basic configuration 
of the ultrasonic diagnosis system. 

Fig. 27 is a diagram which shows a specific 
configuration of an ultrasonic scanner employed in the 
10 ultrasonic diagnosis system. 

Best Mode for Carrying Out the Invention 
Description will be made below regarding an ultrasonic 
diagnosis system according to an embodiment of the present 

15 invention with reference to the accompanying drawings. An 
ultrasonic diagnosis system according to the present 
embodiment employs a method which is referred to as 
"expanded combined autocorrelation method" wherein the 
processing for one -dimensional detection by correlation 

20 calculation based upon the envelope signals using the 

combined autocorrelation method through a one -dimensional 
window is expanded to handle the displacement in the 
horizontal direction by two-dimensional detection through a 
two-dimensional correlation window. Furthermore, with the 

25 expanded combined autocorrelation method according to the 
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present embodiment, envelope -correlation calculation is 
performed for only the grid points set with a pitch of half 
the wavelength of the ultrasonic wave in the axial direction, 
and with the beam- line pitch in the horizontal direction for 
5 reduction of calculation amount, thereby enabling high-speed 
calculation. Note that the expanded combined 
autocorrelation method according to the present embodiment 
employs the phase information in the Scime way as with the 
combined autocorrelation method for improving the precision 

10 of the estimated displacement in the axial direction. 
However, the phase information is not employed for 
estimating displacement in the horizontal direction due to 
lack of the signals serving as a carrier. Accordingly, the 
precision of the estimated displacement in the horizontal 

15 direction is limited corresponding to the sampling pitch 
(ultrasonic -wave beam- line pitch) as with the spatial 
correlation method. Note that the expanded combined 
autocorrelation method according to the present embodiment 
has no particular mechanism for improving the precision of 

20 the estimated displacement in the horizontal direction since 
the elastic modulus distribution can be estimated based upon 
the strain (displacement) distribution in the axial 
direction alone using the elastic modulus distribution 
reconstructing method described later. 

25 Before description of a specific configuration of the 
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expanded combined autocorrelation method according to the 
present embodiment, description will be made regarding the 
combined autocorrelation method which is the basis of the 
expanded combined autocorrelation method according to the 
5 present invention with reference to Fig. 1 through Fig. 6. 
Fig. 1 is a diagram for describing a mechanism of an 
ultrasonic diagnosis apparatus. As can be clearly 
understood from the drawing, an ultrasonic probe 10 serving 
as an ultrasonic detector has functions for converting 

10 electric signals to an ultrasonic wave and converting an 

ultrasonic wave into electric signals, which allows the user 
to cast ultrasonic pulse signals toward the tissue 11. A 
part of the ultrasonic pulse signals passing through the 
tissue 11 is reflected at a first boundary 12 between the 

15 regions having acoustic impedance different one from another. 
The reflected ultrasonic wave which will be referred to as 
"reflected echo signals 12a" passes through the tissue 
toward the ultrasonic probe 10, and the other ultrasonic 
wave passes through the first boundary 12. A part of the 

20 ultrasonic pulse signals passing through the first boundary 
12 is reflected at a second boundary 13 between the regions 
having acoustic impedance different one from another. In 
the same way, the ultrasonic pulse signals reflected at the 
second boundary 13 which will be referred to as "reflected 

25 echo signals 13a" passes through the tissue toward the 
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ultrasonic probe 10, and on the other hand, the other 
ultrasonic pulse signals passes through the second boundary 
13. The reflected ultrasonic echo signals are received by 
the ultrasonic probe 10 so as to be converted into reflected 
5 echo signals which are electric signals. In this case, the 
period of time t from casting of the ultrasonic pulse 
signals from the ultrasonic probe 10 up to reception of the 
echo signals reflected from a reflecting object 14 (boundary 
between the regions having acoustic impedance different one 
10 from another) positioned away from the ultrasonic probe 10 

by the distance L is represented by the following Expression 
(1) . 

<^ ...(1) 

Here, c represents the sound speed within the tissue, 
15 and can be determined to be a constant at around 1500 

[m/ second] through soft tissue. Accordingly, the distance L 
between the probe and the reflecting object is calculated 
based upon the time t from casting of the ultrasonic wave up 
to reception of the reflected echo signals. Furthermore, 
20 the reflected echo signals have information with regard to 
the acoustic properties of the tissue, and accordingly, 
images of the tissue information such as B-mode tomographic 
images can be displayed on a monitor based upon the 
reflected echo signals . 
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For example, a method is known wherein the elastic 
properties which represent the rigidity of the tissue are 
measured using an ultrasonic diagnosis apparatus. The 
aforementioned method for measuring the elastic properties 
5 has a mechanism wherein mechanical vibration is applied to 
the tissue, and the rigidity information is estimated based 
upon the propagating speed of the transverse wave thus 
generated, based upon the fact that the transverse wave 
passes through rigid tissue at a high propagating speed and 

10 passes through soft tissue at a low propagating speed. 

Strictly, the propagating speed v of the transverse wave 
which propagates through the tissue is dependent upon the 
density of the tissue p, the shear modulus |Xi, shear 
viscosity \i2, and the angular frequency of the vibration co, 

15 as represented by the following Expression (2). 



On the other hand, as shown in Fig. 2, a method has 
been proposed wherein static pressure is applied to the 
tissue, and the elastic properties of the tissue are 
20 estimated based upon the strain distribution thereof under 
the applied pressure. The aforementioned method has been 
developed based upon the fact that such static pressure 
causes small strain within rigid tissue, and causes great 




(2) 



strain within soft tissue. Fig. 2(A) is a diagram which 
shows a specific example of a tissue-elasticity measuring 
method with application of static pressure. Fig. 2(B) is a 
diagram which shows the mechanism of the tissue elasticity 
5 measuring method with application of static pressure. As 
can be clearly understood from the drawings, with the 
aforementioned method, normal measurement of the ultrasonic 
echo signals (RF signals without application of pressure) is 
made for the tissue 11 without application of pressure using 

10 a conventional ultrasonic diagnosis apparatus with the 

ultrasonic probe 10. Subsequently, the tissue 11 is pressed 
to a slight degree (around several percent) through the 
ultrasonic probe 10, and measurement of the ultrasonic echo 
signals (RF signals with application of pressure) is made 

15 for the tissue 11 under the applied pressure. Then, the 

displacement distribution which represents the displacement 
of each point within the tissue under pressure is estimated 
based upon the RF signals measured with and without pressure 
applied to the tissue. The principal examples of the 

20 displacement -distribution estimating methods include a 

method using spatial correlation, and a method using the 
Doppler effect. 

Fig. 3 is a diagram which shows a mechanism of the 
spatial correlation method. With the present method, the 

25 strain distribution within the tissue under the applied 
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pressure is estimated based upon the RF signals (or envelope 
signals of the RF signals) with and without pressure applied 
to the tissue by template matching using a two-dimensional 
correlation function. Description will be made below 
5 regarding specific processing thereof. First, with the RF 
signals (or envelope signals thereof) with and without 
pressure applied to the tissue as ii{t, x) and i2(t, x) , 
respectively, the cross -correlation coefficient C(t, x; n, 
m) between the aforementioned two signals is represented by 
10 the following Expression (3). 

/o/2 jfo/2 



I fo/2 V2 I lo/2 Xo/2 



.(3) 



Here, t represents the coordinate point in the 
ultrasonic -wave beam direction (axial direction), x 
represents the coordinate point orthogonal thereto (in the 

15 horizontal direction), tO represents the correlation-window 
size in the axial direction, xO represents the correlation- 
window size in the horizontal direction, Lt represents the 
sampling pitch in the axial direction, and Lx represents the 
sampling pitch in the horizontal direction. Furthermore, n 

20 and m represent integers. In this case, with the 

combination (n, m) of the cross -correlation coefficient 
which exhibits the maximum value as (k, 1) , the displacement 
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Uy in the axial direction and the displacement Ux in the 
horizontal direction at the measurement point (t, x) are 
represented by the following Expressions, respectively. 
Uy = kLt 
5 Ux = ILx 

Note that data sampling is made at a rougher sampling 
pitch Lx in the horizontal direction than the sampling pitch 
Lt in the axial direction, leading to poorer precision of 
the displacement detection in the horizontal direction than 

10 in the axial direction. The aforementioned processing is 
made for each measurement point, whereby displacement 
distribution is estimated. The spatial correlation method 
has the advantage of enabling the user to estimate two- 
dimensional displacement vector components. Furthermore, 

15 the spatial correlation method allows the user to estimate 
the displacement distribution even if great strain (around 
5%) occurs in the tissue. However, the aforementioned 
method leads to a great amount of calculation, leading to 
difficulty in calculation in real time, unlike the 

20 conventional ultrasonic measurement systems. Furthermore, 
the precision of estimated displacement is limited by the 
sampling pitch, thereby leading to a problem of relatively 
poor precision as compared with the Doppler method which 
will be described later. 

25 Fig. 4 is a diagram which shows a mechanism of the 
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Doppler method. With the present method, the displacement 
distribution within the tissue under the applied pressure Is 
estimated based upon the RF signals with and without 
pressure applied to the tissue using the Doppler effect 
5 which Is also employed In blood flow measurement. 

Description will be made below regarding specific processing - 
First, the RF signals with and without pressure applied to 
the tissue are represented by models as represented by the 
following Expression (4). 



Here, li{t) represents the RF signals without 
application of pressure, l2(t) represents the RF signals 
with application of pressure, A{t) represents the envelope 
signals, (Oo represents the center angular frequency of the 
15 ultrasonic wave, and x represents the time shift. Upon 

performing quadrature detection for each of two RF signals, 
the base-band signals are obtained as represented by the 
following Expression. 



10 



...(4) 



'0 



(5) 



20 



Also, the complex correlation function Ri2(t) between 



the aforementioned two signals Is represented by the 
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following Expression. 

£'>,(/+vK(f+v)-rf. = 

Here, RA{t) represents the autocorrelation function of 
the envelope signals, and to represents the correlation 
5 window size in the ultrasonic-beam axial direction. 

Furthermore, the asterisk "*" represents a complex conjugate 
operator. Accordingly, the time shift x and the 
displacement Uy in the axial direction due to application of 
pressure are obtained from the phase (|)(t) of the correlation 
10 function Ri2(t) as represented by the following Expression 
(7) . 

CT 

2 ...(7) 

Note that c represents the sound speed within the 
tissue, and is assumed to be constant within the tissue. 

15 With the Doppler method, the aforementioned processing 

is perfoarmed for each measurement point, and the 
displacement distribution is estimated in the same way as 
with the blood flow measurement developed based upon the 
Doppler effect. Thus, the Doppler method has the advantage 

20 of real-time measurement. Furthermore, the Doppler method 
makes calculation using the phase information, thereby 
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estimating the displacement with higher precision than with 
the spatial correlation method. However, the Doppler method 
has the disadvantage that measurement of large displacement, 
e.g., displacement of a quarter or more the wavelength of 
5 the ultrasonic -wave center frequency leads to aliasing, 

resulting in a problem that correct displacement cannot be 
estimated. Furthermore, the Doppler method has the 
disadvantage that displacement other than that in the axial 
direction cannot be estimated as can be understood from the 

10 above Expression. 

In order to solve the aforementioned problems, the 
present inventors have proposed "Combined Autocorrelation 
Method (CA method)" having the advantages of both of the 
aforementioned two methods. Fig. 5 is a diagram which shows 

15 a mechanism of the combined autocorrelation method proposed 
by the present inventors. With the combined autocorrelation 
method, calculation is made using correlation of the 
envelope signals of the RF signals, thereby solving a 
problem of aliasing which is the disadvantage in the Doppler 

20 method. Description will be made below regarding specific 
processing. 

First, the RF signals with and without pressure applied 
to the tissue are represented by models as represented by 
the following Expression in the same way as with the Doppler 
25 method. 
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A(0 = Re[^(Oe-^^'^"^^ J 

Here, ii(t) represents the RF signals without 
application of pressure, l2(t) represents the RF signals 
with application of pressure, A{t) represents the envelope 
signals, (Oq represents the center angular frequency of the 
ultrasonic wave, and t represents the time shift. Upon 
performing quadrature detection for each of two RF signals, 
the base-band signals are obtained as represented by the 
following Expression. 

Then, the complex correlation function Ri2(t; n) between 
the aforementioned two signals Is represented by the 
following Expression. 

* /2 T ' T -J'^i^-"^ 

(n = .-,-2 "l,OJl,2,"0 .„ ( 10 ) 

15 Here, T represents the cycle of the ultrasonic wave, 

RA(t; t) represents the autocorrelation function of the 
envelope signals, and to represents the correlation window 
size. Furthermore, the asterisk "*" represents a complex 
conjugate operator. Here, n represents the variable number, 

20 and each calculation Is made for different n; displacement 



10 
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at each measurement point being calculated around the point 
determined by the variation number. Note that in a case of 
n = 0, the combined autocorrelation function matches the 
autocorrelation function in the Doppler method as 
5 represented by the Expression (6). That is to say, the 
combined autocorrelation function with n of 0 leads to a 
problem of aliasing in a case wherein measurement is made 
for displacement of a quarter or more the wavelength of the 
ultrasonic wave. In order to solve the aforementioned 
10 problem, with the combined autocorrelation method, the 
envelope correlation coefficient C(t; n) is defined as 
represented by the following Expression (11). 

...(11) 



C7(/;w) = — 7== , , p 

V|A,a;0)|-|J«^C^;/2)l 



Note that Rii(t; 0) represents the autocorrelation 
15 function of Si(t), and R22(t; n) represents the 

autocorrelation function of S2(t + nT/2). With n of the 
envelope correlation coefficient C(t; n) which exhibits the 
maximum value as k, calculation is made using the phase <|) of 
Ri2(t; k) with n k. In this case, displacement is 
20 calculated without aliasing. The reason is that calculation 
of the envelope correlation is made at a pitch of half the 
wavelength. Note that the calculation pitch of half the 
wavelength is the maximum pitch for calculating displacement 
while preventing aliasing. Thus, the time shift x and the 



displacement Uy in the axial direction are calculated using 
<|)(t; k) as represented by the following Expression. 



Note that c represents the sound speed within the 
tissue, and is assumed to be constant within the tissue. 

With the combined autocorrelation method, the 
aforementioned processing is performed for each measurement 
point, whereby displacement distribution is estimated, which 
is an expanded method of the Doppler method. Thus, the 
combined autocorrelation method has the advantage of real- 
time measurement. Furthermore, the combined autocorrelation 
method has the advantage of enabling the user to estimate 
the displacement distribution containing large displacement 
(i.e., displacement of a quarter or more the wavelength of 
the ultrasonic wave) using envelope correlation, unlike the 
Doppler method. 

Fig. 6 is a block diagram which shows a circuit 
configuration for performing basic algorithm of the combined 
autocorrelation method. In the drawing, echo signals x(t) 
obtained without application of pressure are input to a 
unpressed quadrature detection circuit (QD) 131 for 
quadrature detection, and the quadrature -detected signals 




2 



...(12) 
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Ix(t) and Qx(t) thus detected are input to a first 
correlation computing circuit 133 and first correlation 
coefficient computing circuits 1350 through 135N. On the 
other hand, echo signals y{t) obtained under the applied 
5 pressure are input to a first pressed quadrature detection 
circuit (QD) 1320 for quadrature detection, and the 
quadrature -detected signals Y(t) = ly + jQy{Iy(t), Qy(t)) 
thus detected are input to a first correlation coefficient 
computing circuit 1340 and the second correlation computing 

10 circuit 1350. A first delay circuit 134 delays the echo 

signals y(t) by the cycle period T of the ultrasonic wave, 
and the delayed echo signals yl = y(t - T) are input to a 
second pressed quadrature detection (QD) circuit 1321. In 
the same way, a second delay circuit 135 delays the echo 

15 signals yl = y(t - T) which have been delayed by the first 
delay circuit 134 by the cycle period T of the ultrasonic 
wave, and the delayed echo signals y2 = y(t - 2T) are input 
to a next second pressed quadrature detection (QD) circuit 
1322 (not shown). Note that the circuit has N delay 

20 circuits, and the echo signals are consecutively delayed, 
and the echo signals delayed by a value obtained by 
multiplying the cycle period T by an integer are input to 
the corresponding pressed quadrature detection circuit in 
the same way. 

25 The first correlation computing circuit 133 computes a 
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correlation value Rxx based upon the signals Ix and Qx, and 
outputs the correlation value Rxx to second correlation 
coefficient computing circuits 1380 through 138N, The 
second correlation computing circuit 1340 receives the 
5 quadrature -detected signals Iy(t) and Qy(t) from the pressed 
quadrature detection circuit 1320, computes a correlation 
value Ryy based upon the signals ly and Qy, and outputs the 
correlation value Ryy to the second correlation computing 
circuit 1380. The first correlation coefficient computing 

10 circuit 1350 receives the quadrature -detected signals Ix(t) 
and Qx(t) from the unpressed quadrature detection circuit 
131, and the quadrature-detected signals Iy(t) and Qy(t) 
from the first pressed quadrature detection circuit 1320, 
calculates the complex base-band signals Sr and Si based 

15 thereupon, and outputs the base-band signals Sr and Si to a 
third correlation computing circuit 1360 and a phase- 
difference computing circuit 1370. The third correlation 
computing circuit 1360 receives the complex base -band 
signals Sr and Si from the first correlation coefficient 

20 computing circuit 1350, calculates the correlation value 
|Rxy| based thereupon, and outputs the calculated 
correlation value |Rxy| to the second correlation 
coefficient computing circuit 1380. The phase-difference 
computing circuit 1370 receives the complex base -band 

25 signals Sr and Si from the first correlation coefficient 
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computing circuit 1350, and calculates the phase difference 
<J)o(t) based thereupon. The second correlation coefficient 
computing circuit 1380 receives the correlation value Rxx 
from the first correlation computing circuit 133, the 
5 correlation value |Rxy| from the third correlation computing 
circuit 1360, and the correlation value Ryy from the second 
correlation computing circuit 1340, computes the correlation 
coefficient Co(t) based upon the aforementioned correlation 
values, and outputs the calculated correlation coefficient 
10 Co(t). 

The second pressed quadrature detection circuit (QD) 
1321 receives the echo signals Yi = y(t - T) delayed by the 
first delay circuit 134, and outputs the quadrature -detected 
signals Yi(t) = lyi + jQYi (lyi(t), Qyi(t)) to a first 

15 correlation coefficient computing circuit 1341 and a second 
correlation computing circuit 1351. The second correlation 
computing circuit 1341 receives the quadrature-detected 
signals Iyi{t) and Qyi(t) from the second pressed quadrature 
detection circuit (QD) 1321, computes a correlation value 

20 RyiYi based upon the signals lyi(t) and Qyi(t), and outputs 

the correlation value RyiYi to a second correlation computing 
circuit 1381. The first correlation coefficient computing 
circuit 1351 receives the quadrature -detected signals Ix(t) 
and Qx(t) from the unpressed quadrature detection circuit 

25 131, and the quadrature -detected signals lyi(t) and Qyi(t) 
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from the second pressed quadrature detection circuit (QD) 
1321, calculates the complex base-band signals Sri and Sn 
based thereupon, and outputs the base-band signals Sri and 
Sii to a third correlation computing circuit 1361 and a 
5 phase-difference computing circuit 1371. The third 

correlation computing circuit 1361 receives the complex 
base-band signals Sri and Sn from the first correlation 
coefficient computing circuit 1351, calculates the 
correlation value |Rxyx| based thereupon, and outputs the 

10 calculated correlation value |Rxyi| to the second 

correlation coefficient computing circuit 1381. The phase- 
difference computing circuit 1371 receives the complex base- 
band signals Sri and Sn from the first correlation 
coefficient computing circuit 1351, and calculates the phase 

15 difference (t)i(t) based thereupon. The second correlation 

coefficient computing circuit 1381 receives the correlation 
value Rxx from the first correlation computing circuit 133, 
the correlation value |Rxyi| from the third correlation 
computing circuit 1361, and the correlation value RyiYi from 

20 the second correlation computing circuit 1341, computes the 
correlation coefficient Ci(t) based upon the aforementioned 
correlation values, and outputs the calculated correlation 
coefficient Ci(t). 

In the same way, each of the second pressed quadrature- 

25 detection circuit (QD) 1322 through 132N which receive the 
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signals from the corresponding delay circuit downstream from 
the first delay circuit 135, each of the second correlation 
computing circuits 1342 through 134N, each of the first 
correlation coefficient computing circuits 1352 through 135N, 
5 each of the third correlation circuits 1362 through 136N, 

each of the phase-dlf f erence computing circuits 1372 through 
137N, and each of the second correlation coefficient 
computing circuits 1382 through 138N, perform the same 
processing as with the first-stage and second-stage circuit 

10 components as described above, whereby the correlation 

coefficients C2(t) through CN(t), and the phase values <|)2(t) 
through c|)N(t) are output. As described above, the circuit 
for performing the basic algorithm of the combined 
autocorrelation method shown In Fig. 6 has a configuration 

15 wherein the echo signals y(t) under the applied pressure are 
delayed by a period of half the cycle period T/2 of the 
ultrasonic wave (half the wavelength) for each of the delay 
circuits 134 through 13N, and each of the echo signals thus 
delayed are quadrature -detected by the corresponding 

20 quadrature -detection circuit (QD) 1320 through 132N. 

As described above, the strain distribution Is obtained 
by spatial differentiation of the estimated displacement 
distribution under the pressure applied to the tissue. The 
strain distribution represents the relative elastic property 

25 of the tissue, and accordingly, diagnosis based upon the 
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strain distribution exhibits effects similar to those of the 
diagnosis based upon the elastic modulus distribution. 
However, in a case of liver cirrhosis which causes rigidity 
of the entire affected tissue, it is difficult to make the 
5 same diagnosis of the tissue as with the elastic-property 
distribution which allows the surgeon to make quantitative 
estimation. Accordingly, in recent years, reconstruction 
methods for tissue elastic modulus distribution are being 
studied. Note that all of these reconstruction methods are 

10 in the research stage, and that no standard method has been 
established. 

On the other hand, the tissue elastic modulus 
distribution can be obtained based upon the strain 
distribution and the stress distribution within the tissue 

15 as described above. However, it is difficult to make direct 
measurement of the stress distribution with the existing 
technique. Accordingly, in practicality, the elastic 
modulus distribution is reconstructed based upon the strain 
distribution so as to satisfy the boundary conditions under 

20 the pressure applied to the tissue, i.e., there is the need 
to solve the inverse problem. In general, it is difficult 
to solve the inverse problem, and few elastic modulus 
reconstruction methods have been proposed. Description will 
be made regarding the conventional elastic modulus 

25 reconstruction methods. 
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First, a method is known, which has been proposed on 
the assvimption that the tissue is represented by a one- 
dimensional model (one -dimensional elastic model). That is 
to say, a method is known wherein the elastic modulus 
5 distribution is calculated on the assumption that the tissue 
is represented by a one -dimensional elastic model. On the 
aforementioned assumption, the elastic modulus is determined 
to be the inverse number of the strain. Strictly, the 
aforementioned method is not an elastic modulus 

10 reconstruction method. With the aforementioned method, only 
the inverse number of strain is obtained, and accordingly, 
only relative elastic property of the tissue can be obtained 
as with the strain distribution. 

Second, a method is known wherein elasticity equation 

15 is reduced to that without the stress terms (on the 

assumption that the tissue exhibits isotropic elasticity, 
incompressibility , and plane strain). With the 
aforementioned method, the elasticity equation formed so as 
to represent the plane- stress state is reduced to that 

20 without the stress terms, and the tissue elastic modulus 
distribution is reconstructed based upon the strain 
distribution (all the components of the strain tensor 
including the shear strain components) using the reduced 
elasticity equation without the stress terms so as to 

25 satisfy the boundary conditions (applied-pressure 
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distribution on the body surface, or the displacement 
distribution thereof). Note that the present method 
requires the region (reference region) where the elastic 
modulus has been obtained beforehand. 
5 Third, a method is known wherein the elasticity 

differential equation is integrated (on the assumption that 
the tissue exhibits isotropic elasticity, incompressibility , 
and plane strain). With the aforementioned method, the 
elasticity equation formed so as to represent the plane- 

10 stress state is reduced to that without the stress terms, 

and the tissue elastic modulus distribution is reconstructed 
based upon the strain distribution (all the components of 
the strain tensor including the shear strain components) by 
consecutively integrating the reduced differential equation 

15 without the stress terms with regard to the elastic modulus 
with the elastic modulus near the body surface as the 
reference. Note that the present method requires the region 
(reference region) near the body surface where the elastic 
modulus distribution has been obtained beforehand. 

20 Furthermore, the aforementioned method has a problem that 
the error of calculation becomes greater the farther away 
from the body surface, due to accumulation of the error 
since integration is made with the elastic modulus near the 
body surface as the reference. 

25 Fourth, a method using the perturbation method is known 
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(on the assumption that the tissue exhibits isotropic 
elasticity, near-incompressibility , and plane strain) . With 
the aforementioned method, the tissue elastic modulus 
distribution is reconstructed by solving the elasticity 
5 equation which has been formed so as to represent the plane- 
strain state, based upon the applied-pressure distribution 
on the body surface and the strain distribution thereof in 
the ultrasonic-beam direction (axial direction) using the 
perturbation method with the iterative method. 

10 Description has been made regarding the basic mechanism 

and the specific problems which are to be solved. 
Description will be made below regarding an embodiment 
according to the present invention, in order to solve the 
aforementioned problems. Fig. 7 is a block diagram which 

15 shows a schematic configuration of an ultrasonic diagnosis 
system according to an embodiment of the present invention. 
An ultrasonic probe 91 comprises a conventional sector scan 
probe (sector phased array probe), a linear scan probe 
(linear array probe) , or a convex scan probe (convex array 

20 probe), having functions for casting an ultrasonic wave 

toward the tissue which is to be observed, and receiving the 
reflected ultrasonic wave. 

The RF signals obtained with and without pressure 
applied to the tissue are output from the ultrasonic probe 

25 91 to a quadrature detector 92. The quadrature detector 92 



- 39 - 



converts the RF signals with and without pressure applied to 
the tissue into the complex envelope signals (IQ signals) 
with and without application of pressure, and outputs the IQ 
signals to a complex two-dimensional correlation calculation 
5 unit 93- The complex two-dimensional correlation 

calculation unit 93 calculates two-dimensional correlation 
between the RF signals with and without pressure applied to 
the tissue, outputs the position which exhibits the maximum 
correlation to a horizontal-direction displacement 

10 calculation unit 94 and an axial -direct ion displacement 

calculation unit 95, and outputs the corresponding phase of 
the correlation function to an axial -direct ion calculation 
unit 95. Note that correlation is calculated in the axial 
direction at a pitch of a half the ultrasonic center 

15 frequency which is the maximum pitch for obtaining the phase 
while preventing aliasing. Correlation is calculated at 
such a pitch for enabling real-time display of the 
ultrasonic diagnosis system. Accordingly, the present 
invention is not restricted to calculation at a pitch of a 

20 half the wavelength, rather, an arrangement may be made 
wherein correlation is calculated with high precision. 

The horizontal-direction displacement calculation unit 
94 calculates the displacement Ux in the horizontal 
direction based upon the position corresponding to the 

25 maximum correlation in the horizontal direction received 
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from the complex two-dimensional correlation calculation 
unit 93, and outputs the displacement to a horizontal- 
direction strain calculation unit 96. On the other hand, 
the axial-direction displacement calculation unit 95 
5 calculates the displacement Uy in the axial direction based 
upon the position corresponding to the maximum correlation 
in the axial direction and the phase received from the 
complex two-dimensional correlation calculation unit 93, and 
outputs the displacement to an axial -direction strain 

10 calculation unit 97. The horizontal-direction strain 

calculation unit 96 makes spatial differentiation of the 
displacement Ux in the horizontal direction received from 
the horizontal-direction calculation unit 94 so as to 
calculates the strain distribution Ex in the horizontal 

15 direction, and outputs the strain distribution to a 
quantization unit 98. On the other hand, the axial- 
direction strain calculation unit 97 makes spatial 
differentiation of the displacement Uy in the axial 
direction received from the axial -direction calculation unit 

20 95 so as to calculate the strain distribution Ey in the axial 
direction, and outputs the strain distribution to the 
quantization unit 98. The quantization unit 98 quantizes 
the strain distribution Ex in the horizontal direction and 
the strain distribution Ey in the axial direction in order to 

25 make grayscale display (or color display) thereof, and 
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displays the information on a display unit 99. The display 
unit 99 displays each of the strain distributions thus 
quantized • 

Next, description will be made regarding the operation 
5 of the expanded combined autocorrelation method employed in 
the ultrasonic diagnosis system shown in Fig. 7. First, let 
us consider a case wherein the tissue is compressed to a 
slight degree (i.e., several percent or less). In this case, 
from the local perspective, the pressure is assximed to cause 
10 parallel displacement for each point within the tissue. 
That is to say, the RF signals with and without pressure 
applied to the tissue are represented by the model 
represented by the following Expression. 

£,(r,x) =ReU(/-r,x-«,)^-^l-('-'>-''lJ ^3 ^ 

15 Here, ii(t, x) represents the RF signals without 

application of pressure, izCt, x) represents the RF signals 
under pressure, A{t, x) represents the envelope signals, coo 
represents the center angular frequency of the ultrasonic 
wave, T represents the time shift serving as a time 

20 parameter which represents the displacement in the axial 

direction, Ux represents the displacement in the horizontal 
direction, and 0 represents the initial phase. Note that 
with the present method, the RF signals with and without 
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application of pressure are represented by a model giving 
consideration to the displacement in the horizontal 
direction, unlike the Doppler method and the combined 
autocorrelation method. The parameter which is to be 
5 obtained in the final stage is the displacement Uy = cx/2 in 
the axial direction (i.e., the time shift x) and the 
displacement Ux in the horizontal direction. Note that c 
represents the sound speed within the tissue, and is assumed 
to be constant within the tissue. 

10 Then, the RF signals with and without pressure applied 

to the tissue are quadrature-detected by the quadrature 
detector 92. That is to say, the sine wave and the cosine 
wave having the same frequency as the center frequency of 
the ultrasonic wave are applied to the RF signals, following 

15 which low-pass filtering is further applied to the RF 

signals, whereby the complex base-band signals Si and S2 are 
obtained as represented by the following Expression (14). 



s^it,xy = ^(/-r.x-wje^^-^'-"' ^^^^^ 



Then, the two-dimensional complex correlation function 
20 Ri2{t, x; n, m) between the Si(t, x) and S2(t + nT/2, x + mL) 
is defined as represented by the following Expression (15). 
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Si* 

(«= - AT^, -.-2,-1,04,2, -sAU) 

(m=-M^,--,-2,.l,0,l,2,-,M„) ^^^^j 

Here, T represents the cycle of the ultrasonic wave, L 
represents the sampling pitch (beam-line pitch), RA(t, x; x, 
Ux) represents the autocorrelation function of the envelope 
5 signals, tO represents the length of the correlation window 
in the axial direction, xO represents the length of the 
correlation window in the horizontal direction. On the 
other hand, v represents a variable value in the time (x) 
axial direction for integration, and w represents a variable 

10 value in the beam-line direction for integration, and the 

asterisk "*" represents a complex conjugate operator. Then, 
the two-dimensional envelope correlation coefficient C(t, x; 
n, m) is defined using the two-dimensional complex 
correlation function as represented by the following 

15 Expression (16). 

...(16) 



Note that Rii(t, x; 0, 0) represents the autocorrelation 
function of Si(t, x) , and R22(t, x; n, m) represents the 
autocorrelation function of S2(t + nT/2, x + mL). The 



envelope correlation coefficient is used for solving a 
problem of aliasing in the same way as with the combined 
autocorrelation method. That is to say, combinations {C{t, 
x; n, m) , (t)(t, x; n, m)} formed of C(t, x; n, m) and <|){t, x; 
5 n, m) of Ri2(t, x; n, m) are obtained for all the variable 
numbers n and m at each measurement point (t, x) . With the 
present method, let us say that the variable number pair (n, 
m) is determined in a sufficient range, i.e., in a range 
sufficient for envelope correlation. In this case, the 

10 phase (|)(t, x; k, 1) corresponding to (n, m) = (k, 1) which 
exhibits the maximum envelope correlation coefficient 
matches the phase without aliasing. The reason is that with 
the variable number pair (n, m) which exhibits the maximum 
envelope correlation coefficient as (k, 1), the time shift 

15 |t - kT/2| between Si{t, k) and SaCt + kT/2, x + IL) is 

smaller than T/2. That is to say, |(l)(t, x; k, 1) | = (o© 1 1 - 
kT/2 1 is smaller than Jt. That is to say, with the present 
method, calculation is made using <|)(t, x; k, 1) without 
aliasing, thereby obtaining the correct time shift t, the 

20 correct displacement Uy in the axial direction, and the 

correct displacement Ux in the horizontal direction, at each 
measurement point (t, x) as represented by the following 
Expression ( 17 ) - 
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^ <f>(t.x;k,/) ^^r 

c r 

^x^^ ...(17) 

Note that c represents the sound speed within the 
tissue (which is assumed to be constant at 1500 m/s which is 
normal sound speed within soft tissue). With the present 
5 method, the calculation described above is made for all the 
measurement points, thereby obtaining the displacement 
distribution Uy(x, y) in the axial direction and Ux(x, y) in 
the horizontal direction. 

Furthermore, with the present method, spatial 
10 differentiation is made for each of the aforementioned 

displacement distributions, whereby the strain distribution 
Ey(x, y) in the axial direction and the strain distribution 
ex(x, y) in the horizontal direction are obtained as 
represented by the following Expression (18). 



dy 

15 ..-(18) 

As described above, with the present method, the 
displacement (strain) distribution in the axial direction 
and in the horizontal direction is estimated based upon the 
RF signals with and without pressure applied to the tissue. 

20 Note that as can be understood from the above Expression Ux 
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= IL, the precision of the displacement detection in the 
horizontal direction is limited by the sampling pitch (beam- 
line pitch) in the horizontal direction, and accordingly, 
the present method has the disadvantage of relatively low 
5 precision in the horizontal direction. However, the present 
method has the advantage of real-time observation, thereby 
improving practical performance . 

The expanded combined autocorrelation method described 
above has a function for analyzing the relative displacement 

10 in the horizontal direction between the tissue and the 

ultrasonic probe under the pressure applied to the tissue 
using a two-dimensional correlation window applied in a 
predetermined range for each calculation, thereby handling 
displacement of the tissue in the horizontal direction. 

15 However, the present two-dimensional expanded combined 
autocorrelation method having such a function cannot 
estimate strain distribution containing displacement in the 
direction (direction orthogonal to the two-dimensional 
ultrasonic scanning plane (slice direction)) orthogonal to 

20 both the axial direction and the horizontal direction under 
the pressure applied to the tissue. In order to solve the 
aforementioned problem, the above two-dimensional expanded 
combined correlation method is easily expanded into the 
three-dimensional expanded combined autocorrelation method 

25 using a three-dimensional window applied to a three- 
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dimensional range, thereby enabling the system to be more 
stable. 

Fig. 9 and Fig. 10 are flowcharts for describing basic 
algorithm of the three-dimensional combined autocorrelation 
5 method. Note that Fig. 10 is a flowchart for making 

detailed description with regard to a part of the processing 
shown in Fig. 9. 

In Step Sll, a scanning-line resistor 1 stores "1", 
which serves as a counter for performing the same processing 
10 for the first scanning line through the L'th scanning line. 
The value stored in the scanning- line resistor 1 is checked 
in determination processing in Step S23. 

In Step S12, the variable in the thickness direction 
(frame direction) is incremented in a range between -U and U 
15 for each loop processing. Note that the count is checked in 
the determination processing in Step S18. 

In Step S13, the variable in the scanning direction 
(scanning line direction) is incremented in a range between 
-V and V for each loop processing. Note that the count is 
20 checked in the determination processing in Step S17. 

In Step S14, the variable in the distance direction 
(axial direction) is incremented in a range between 0 and M 
for each loop processing. Note that the count is checked in 
the determination processing in Step S16. 
25 In Step S15, correlation coefficient C(l, t; u, v, n) 
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of the envelope signals in the distance direction (axial 
direction) is calculated with the combined autocorrelation 
method. Note that the conventional combined autocorrelation 
method cannot exhibit sufficient calculation speed for 
5 three-dimensional calculation, and accordingly, the 

correlation coefficient C(l, t; u, v, n) is calculated with 
the high-speed combined autocorrelation method. Description 
will be made later regarding the high-speed autocorrelation 
method. 

10 In Step S16, processing is performed for the variable 

determined in the aforementioned Step S14. That is to say, 
determination is made whether or not the variable n stored 
in the distance-direction register has reached the 
predetermined maximum M. In the event that determination is 

15 made that the variable n has reached M, the flow proceeds to 
Step S17. Otherwise, the flow returns to Step S14, and the 
variable n stored in the distance-direction register is 
incremented . 

In Step S17, processing is performed for the variable 
20 determined in the aforementioned Step S13. That is to say, 
determination is made whether or not the variable v stored 
in the scanning- direction register has reached the 
predetermined maximum V. In the event that determination is 
made that the variable v has reached V, the flow proceeds to 
25 Step S18. Otherwise, the flow returns to Step S13, and the 
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variable v stored in the scanning-direction register is 
incremented . 

In Step S18, processing is performed for the variable 
determined in the aforementioned Step S12. That is to say, 
5 determination is made whether or not the variable u stored 
in the thickness-direction register has reached the 
predetermined maximiim U. In the event that determination is 
made that the variable u has reached U, the flow proceeds to 
Step S19. Otherwise, the flow returns to Step S12, and the 
10 variable u stored in the thickness-direction register is 
incremented . 

In Step S19, the system searches the variable 
combination (u, v, n) in a range of u = {-U, 0, U) , v 

= (-V, 0, V), and n = (0, 1 N) for the variable 

15 combination (uo, vo, n©) which exhibits the maximum 

correlation coefficient C(l, t; u, v, n) calculated in Step 
S12 through Step S18. Note that while the present 
embodiment employs such a distance-direction range n = (0, 1, 
N) under the assumption that only the positive pressure 
20 is applied to the tissue, it is needless to say that the 
distance-direction range n = (-M, 0, N) should be 
employed in an arrangement wherein both the positive and 
negative pressure is applied to the tissue. 

In Step S20, the phase difference (|)(1, t; Uo, Vo, no) of 
25 the correlation coefficient C(l, t; Uo, Vo, no) obtained in 
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Step S19 is calculated. 

In Step S21, the phase difference noJt + (|)(1, t; Uo, Vo, 
no) is calculated as the phase difference in the final stage. 

In Step S22 , the displacement v = Vo + Av in the 
5 scanning direction, and u = Uo + Au in the thickness 

direction, are calculated using the correlation coefficient 
C(l, t; u, V, n) at a point (u, v, n) near the point (uo, v©, 
no) . 

In Step S23, the variable stored in the scanning-line 
10 resistor 1 in the aforementioned Step Sll is checked. That 
is to say, determination is made whether or not the variable 
1 has reached L. In the event that determination is made 
that 1 has reached L, the flow proceeds to Step S24. 
Otherwise, the flow returns to Step Sll, and the variable 
15 stored in the scanning- line resistor 1 is incremented. 

In Step S24, strain distribution is calculated by 
making spatial differentiation of the estimated displacement 
distribution under pressure applied to the tissue. 

Fig. 10 is a flowchart for making detailed description 
20 regarding the aforementioned high-speed combined 

autocorrelation method in Step S15 shown in Fig. 9. 

In Step S31, the RF signals x without application of 
pressure and the RF signals y under pressure are quadrature- 
detected, whereby I signals and Q signals are obtained as 
25 described below. 
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x(t) -> Ix, Qx (where X(t) = Ix + jQx) 
y{t) ly, Qy (where Y(t) = ly + jQy) 
In Step S3 2, the correlations Rxy, Rxx, and Ryy, are 
calculated based upon the following Expressions. 
5 Rxy = / X(t + V) • Y*(t + v) dv 

Rxx = / X(t + V) • X*(t + V) dv 
Ryy = / Y(t + V) • Y*(t + V) dv 

In Step 833, the correlation coefficient Co is 
calculated based upon the correlations Rxy, Rxx, and Ryy, 
10 thus obtained, using the following Expression. 
Co = |Rxy|/(VRxx ' VRyy) 

In Step S34, the variable number n is set to 1. 
In Step S35, Yn(t) = Y(t - nT) e^'"'''^ is calculated. 
In Step S36, Rxyn and RynYn are calculated based upon 
15 the following Expressions. 

Rxyn = / X(t + V) • Yn*(t + V) dv 

= / X(t + V) • Y*(t - nT + V) e^'^'''^ dv 

RynYn - / Yn (t + V) • Yn* ( t + V) dV 

= / Y(t -nT+v) • Y*(t-nT+v) dv 
20 In Step S37, the correlation coefficient Cn is 

calculated based upon the Rxyn and RynYn thus obtained, using 
the following Expression. 

Cn = |Rxyn| /(v^Rxx • VRynYn) 

In step S38, the variable number n is incremented. 
25 In step S39, determination is made whether or not the 
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variable number n has reached the predetermined maximum 
number M. In the event that determination is made that n 
has reached M, the processing ends. Otherwise, the flow 
returns to Step S35, and the same processing is repeated. 
5 With the method shown in the flowchart in Fig. 10, Yn is 

calculated from Y in Step S35 in order to calculate Rxyn and 
RYnYn- This enables a simple circuit configuration. 
Description will be made below regarding a method for 
calculating Yn from Y. 
10 First, the echo signals x(t) without application of 

pressure are represented by the following Expression. 
x(t) = u(t) cos((ot + 6) 

On the other hand, the echo signals y(t) under pressure 
applied in the cLxial direction are represented by the 
15 following Expression. 

y(t) = x(t + T) = u(t + x) cos(co(t + x) + 6). 

Then, the quadrature -detected signals of x, y, and yn, 
are represented by the following Expressions. 

x(t) ^ 

20 Ix = 0.5 u(t) cosO 

Qx = -0.5 u(t) sine 
(X{t) = Ix + jQx = 0.5 u(t) e-^^) 

y(t) ^ 

ly = 0.5 u(t + x) cos((ox + 9) 
25 Qy = -0.5 u(t + x) sin(cax + 9) 
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(Y(t) = ly + jQy = 0.5 u(t + t) e"^^^ * ®M 
yn(t) = y{t - nT) 

= u(t + X - nT) cos(co{t + T - nT) + 6) 
= u(t + T - nT) cos(a)t + (o(t - nT) + 9) 
5 Here, T represents the half cycle. Then, lyn and Qyn 

are obtained as represented by the following Expressions, 
lyn = 0.5 u(t + T - nT) cos(co(t - nT) + 0) 
Qyn =-0.5u(t+T-nT) sin{co(T - nT) +9) 
(Yn ^ lyn + JQyn = 0.5 u(t + T - nT ) e-^^"<^ " '^'^> ^ 
10 Accordingly, the following relation is obtained based 

upon the above Expressions. 

yn(t) = lyn + jQyn 

= 0.5 u(t + T - nT) e-^^"«" " ^''^ ^ 
= Y(t - nT) e^"""" 

15 As can be understood from the above Expression, Yn(t) is 

calculated from Y(t) = ly + jQy. 

Accordingly, Rxyn and RynYn are calculated from X and Y 
as represented by the following Expressions. 
Rxyn = 4/ X{t + V) • Yn*{t + V) dv 
20 4/ X(t + V) • Y*{t - nT + V) e^'"'''^ dv 

I Rxyn I = RUn 

= / u(t+v)u(t+x - nT + v) dv 

=4/ |X(t + V) • Yn*(t + v)| dv 

= 4/ |X(t + V) • Y*(t - nT + V) e^'^'''^ | dv 

25 = 4/ |X(t + v) • Y*{t - nT + v) I dv 
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RynYn = / u(t+T-nT + v)u(t+x-nT + v)dv 

= 4/ |Yn{t + V) • Yn*(t + V) I dV 

= 4/ Y(t - nT + v) • Y*(t - nT + v) dv 
Here, the asterisk "*" represents a complex conjugate 
5 operator . 

Fig. 11 is a block diagram which shows a circuit 
configuration for executing the basic algorithm of the 
three-dimensional combined autocorrelation method. With the 
conventional circuit configuration for executing the 

10 combined autocorrelation method as shown in Fig. 7, a great 
number of quadrature-detection circuits 1320 through 132N 
are arrayed, and processing in such a number of the 
quadrature-detection circuits 1320 through 132N requires 
enormous processing time, leading to difficulty in high- 

15 speed calculation, resulting in difficulty in real-time 

image display. Accordingly, the present embodiment employs 
the circuit configuration for executing the above -described 
basic algorithm as shown in Fig. 11, thereby enabling a 
high-speed circuit for executing the combined 

20 autocorrelation method. 

The unpressed quadrature- detection circuit (QD) 131 
receives the echo signals x(t) without application of 
pressure, quadrature-detects the received echo signals, and 
outputs the quadrature -detected signals Ix(t) and Qx(t) to 

25 the first correlation computing circuit 133 and the first 
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correlation coefficient computing circuits 1350 through 135N. 
The pressed quadrature-detection circuit (QD) 132 receives 
the echo signals y(t) under pressure, quadrature -detects the 
received echo signals, and outputs the quadrature -detected 
5 signals Y(t) = ly + jQy (Iy(t), Qy(t)) to the first 

correlation coefficient computing circuits 1350, the second 
correlation computing circuit 1340, the first delay circuit 
134, and second delay circuit 135. The first delay circuit 
134 and the second delay circuit 135 delay the quadrature - 

10 detected signals Y(t) by the cycle T of the ultrasonic wave, 
and output the delayed quadrature -detected signals Y(t - T) 
to the first correlation computing circuit 1351, the third 
delay circuit 136, and the fourth delay circuit 137. Each 
of the third delay circuit 136 and the fourth delay circuit 

15 137 delay the quadrature -detected signals Y(t - T) by the 
cycle T of the ultrasonic wave, and output the delayed 
quadrature -detected signals Y(t - 2T) to the subsequent 
first correlation coefficient computing circuit and delay 
circuit (not shown) . In the same way, the quadrature - 

20 detected signals are consecutively delayed by the cycle T 
using each of the N delay circuit, and the delayed signals 
are output to the corresponding first correlation 
coefficient computing circuit. 

The first correlation computing circuit 133 computes 

25 the correlation value Rxx based upon the signals Ix and Qx, 
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and outputs the correlation value Rxx to the second 
correlation coefficient computing circuits 1380 through 138N. 
The second correlation computing circuit 1340 receives the 
quadrature -detected signals Iy(t) and Qy(t) from the pressed 
5 quadrature -detect ion circuit 132, computes the correlation 
value Ryy based upon the signals ly and Qy, and outputs the 
correlation value Ryy to the second correlation coefficient 
computing circuit 1380. The first correlation coefficient 
computing circuit 1350 receives the quadrature-detected 

10 signals Ix(t) and Qx(t) from the unpressed quadrature- 
detection circuit 131, and the quadrature -detected signals 
Iy(t) and Qy(t) from the pressed quadrature-detection 
circuit 132, calculates the complex base-band signals Sr and 
Si, and outputs the complex base-band signals Sr and Sj to 

15 the third correlation computing circuit 1360 and the phase- 
difference computing circuit 1370. The third correlation 
computing circuit 1360 receives the complex base-band 
signals Sr and Si from the first correlation coefficient 
computing circuit 1350, calculates the correlation value 

20 |Rxy| based thereupon, and outputs the correlation value 
|Rxy| to the second correlation coefficient computing 
circuit 1380. The phase-dif f erence computing circuit 1370 
receives the complex base-band signals Sr and Si from the 
first correlation coefficient computing circuit 1350, and 

25 calculates the phase difference (|)o(t) based thereupon. The 
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second correlation computing circuit 1380 receives the 
correlation value Rxx from the first correlation computing 
circuit 133, the correlation value |Rxy| from the third 
correlation computing circuit 1360, and the correlation 
5 value Ryy from the second correlation circuit 1340, computes 
the correlation coefficient Co(t) based upon these 
correlation values, and outputs the correlation coefficient 
Co(t) . 

The second correlation computing circuit 1341 receives 

10 the delayed quadrature -detected signals Iy(t - T) and Qy(t - 
T) from the first delay circuit 134 and the second delay 
circuit 135, computes the correlation value Ryiyi based upon 
the signals Iy(t - T) and Qy(t - T) , and outputs the 
correlation value Ryiyi to the second correlation coefficient 

15 computing circuit 1381. The first correlation coefficient 
computing circuit 1351 receives the quadrature -detected 
signals Ix(t) and Qx(t) from the unpressed quadrature- 
detection circuit 131, and the delayed quadrature -detected 
signals Iy(t - T) , and Qy(t - T) from the first delay 

20 circuit 134 and the second delay circuit 135, obtains the 
complex base -band signals Sri and Sn, and outputs the 
complex base -band signals Sri and Sn to the third 
correlation computing circuit 1361 and the phase-difference 
computing circuit 1371, The third correlation computing 

25 circuit 1361 receives the complex base-band signals Sri and 
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Sii from the first correlation coefficient computing circuit 
1351, obtains the correlation value | Rxyi | based thereupon, 
and outputs the correlation value |Rxyi| to the second 
correlation coefficient computing circuit 1381. The phase- 
5 difference computing circuit 1371 receives the complex base- 
band signals Sri and Sn from the first correlation 
coefficient computing circuit 1351, and obtains the phase 
difference <l)i(t) based thereupon. The second correlation 
coefficient computing circuit 1381 receives the correlation 

10 value Rxx from the first correlation computing circuit 133, 
the correlation value |Rxyi| from the third correlation 
computing circuit 1361, and the correlation value RyiYi from 
the second correlation computing circuit 1341, calculates 
the correlation coefficient Ci(t) based these correlation 

15 values, and outputs the correlation coefficient Ci{t). 

In the same way, with the second correlation computing 
circuits 1342 through 134N, the first correlation 
coefficient computing circuits 1352 through 135N, the third 
correlation computing circuits 1362 through 136N, the phase- 

20 difference computing circuits 1372 through 137N, and the 
second correlation coefficient computing circuits 1382 
through 138N, arrayed downstream from the third delay 
circuit 135 and the fourth delay circuit 136, the same 
processing as described above is performed for each of the 

25 consecutively delayed quadrature -detected signals Iy(t - 2T) , 
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Iy(t - NT), and Qy(t - 2T) , Qy(t - NT), whereby the 
correlation coefficients CaCt) through CwCt), and the phase 
(t)2 { t ) through 4>n ( t ) , are output . 

Next, description will be made regarding the elastic- 
5 modulus distribution reconstructing method using the three- 
dimensional finite element method. In order to simplify the 
inverse problem for reconstructing the elastic modulus 
distribution, with the present embodiment, the tissue is 
represented by a model. This allows the user to perform the 

10 elastic-modulus distribution reconstructing method proposed 
in the present embodiment using the finite element method. 
Specifically, with the present embodiment, the tissue is 
represented by a model on the assumption described below. 
First, the tissue is assumed to exhibit isotropic 

15 elasticity. On the other hand, the distribution of the 
tissue strain is estimated for the tissue under external 
static pressure. Note that the static pressure is applied 
to the tissue so as to cause slight compression thereof, 
thereby allowing the user to make calculation using the 

20 relation between the RF signals with and without pressure 
applied to the tissue. Accordingly, in this case, the 
relation between the stress and the strain exhibits 
linearity. That is to say, approximate calculation can be 
made on the assumption that the tissue is represented by an 

25 elastic model. Furthermore, in this case, the tissue is 
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assumed to be isotropic, and accordingly, the tissue is 
assumed to exhibits isotropic elasticity. 

Furthermore, the tissue is assumed to exhibit near- 
incompressibility . In general, it is known that the tissue 
5 exhibits compressibility near incompressibility (Poisson's 
ratio V = 0.5). With the present embodiment, calculation is 
made with a constant Poisson's ratio of 0.49 within the 
tissue. Note that with the present embodiment, calculation 
is not made on the assumption that the tissue exhibits the 

10 complete incompressibility. The reason is that if 

calculation is made on the asstmiption that the tissue 
exhibits the complete incompressibility, i.e., calculation 
is made with a constant Poisson's ratio of 0.5 within the 
tissue, the elastic equation becomes a special case, leading 

15 to a problem that calculation cannot be made using the 

finite element method proposed in the present embodiment . 
Furthermore, with the present embodiment, the Poisson's 
ratio is assumed to be constant within the tissue, and 
accordingly, only the Young's modulus should be estimated 

20 for the elastic -modulus distribution, thereby reducing the 
inverse problem. Note that in general, the Poisson's ratio 
exhibits a relatively small variation as compared with the 
Young's modulus. Accordingly, with the present embodiment, 
calculation is made with the constant Poisson's ratio of 

25 0.49 within the tissue. 
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Next, the tissue is represented by a three-dimensional 
finite element model. With the reconstructing method for 
the elastic modulus distribution according to the present 
embodiment, calculation is made using the finite element 
5 method, and accordingly, the tissue is divided into a finite 
number of rectangular parallelepiped elements. Note that 
each element is assumed to exhibit the constant elastic 
modulus, the constant stress, and the constant strain, 
therewithin. In general, in order to understand the method 

10 for solving the inverse problem, it is important to 

understand the forward problem corresponding to the inverse 
problem. With the present embodiment, the inverse problem 
is that the elastic modulus distribution is estimated based 
upon the strain distribution and the boundary conditions. 

15 Accordingly, the forward problem corresponding to the 
aforementioned inverse problem is that the strain 
distribution is calculated based upon the elastic modulus 
distribution and the boundary conditions. With the present 
embodiment, the aforementioned problem is solved using the 

20 finite element method (FEM: Finite Element Method) which is 
a kind of known numerical analysis methods. 

With the finite element method, the tissue serving as a 
continuum which is to be estimated is represented by a model 
formed of a finite number of elements, and simultaneous 

25 linear equations which represent the properties of each 
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element are solved using numerical analysis. Note that 
description will be made later regarding specific 
calculation using the finite element method. In brief, with 
the finite element method, the strain (displacement) 
5 distribution and the stress distribution serving as the 
output values are obtained based upon the elastic -modulus 
distribution of the tissue and the boundary conditions 
serving as the input values. 

With the present embodiment, approximate calculation is 

10 made on the assumption that the tissue exhibits isotropic 
elasticity, and accordingly, the elastic equations 
(equilibrium equation, relation between strain and 
displacement, relation between stress and strain) hold under 
the conditions within the tissue as represented by the 

15 following Expressions. 

The equilibrium equation is represented by the 
following Expression (19). 

(/,y = l,2^) ...(19) 

Here, the reference numeral 1, 2, and 3, serving as i 
20 and j, represent x, y, z, respectively. On the other hand, 
the relation between strain and displacement is represented 
by the following Expression (20). 




...(20) 



The relation between stress and strain (generalized 
Hooke's law) is represented by the following Expression (21). 



The aforementioned elastic equations are represented 
using tensors. Accordingly, there are actually three 
equilibrium equations, six strain-displacement relations, 
and six stress -strain relations. Note that the coordinate 
(xl, x2, x3) represents (x, y, z). Other reference 
characters represent the properties as follows. 

E: Young's modulus (which is also referred to as 
"elastic modulus") 

V: Poisson's ratio 

Eij: strain tensor 

(enn = Ell + E22 + £33: strain of volume) 
Si-j: stress tensor 
6ij : Kronecker delta 
Ui: displacement vector 
fi: volume force vector 

(Note that the gravity is deemed negligible, and 
accordingly, fi is assumed to be zero in this case) 

The relation between the stress and the strain is then 
solved for eij. As a result, the relation between strain and 




(21) 
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Stress is transformed as represented by the following 



Expression ( 22 ) . 



...(22) 



Here, Snn = Sn + S22 + 533- Then, i = j = 2 is 

5 substituted into the above Expression (22), and the 

Expression is solved for the Young's modulus E, thereby 
obtaining the following Expression (23). 



10 the Young's modulus, i.e., the elastic modulus can be 

calculated based upon the strain component in the axial 
direction (ultrasonic-wave beam direction, in the present 
embodiment), and the stress components in all directions. 
However, it is difficult to make direct measurement of the 

15 stress distribution used for the above Expression with the 
current techniques. Accordingly, with the present 
embodiment, the stress distribution and the elastic-modulus 
distribution are alternately estimated and updated such that 
the estimated elastic modulus distribution becomes close to 

20 the actual distribution. The specific procedure for 

reconstructing the elastic modulus distribution is performed 




22 



...(23) 



As can be understood from the above Expression (23), 
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as follows . 

First, let us say that a uniform distribution is 
employed as the initial distribution {e"^} for estimating the 
elastic modulus distribution. Second, the stress 
5 distribution {S^> caused due to the initial elastic modulus 
distribution {e'°} is calculated using the three-dimensional 
finite element method. Specifically, the strain- 
displacement relation and the stress -strain relation are 
substituted into the above equilibrium equation, whereby a 
10 new equilibrium equation is formed as represented by the 

following Expression (24). The new equilibrium equation is 
applied to each element within the tissue model. 

Here , 

15 9x„ dx^ dx^ dx^ ...(25) 

Then, the simultaneous equations are solved for the 
displacements using the Gaussian elimination under the 
following boundary conditions, whereby the displacement 
distribution {u corresponding to the elastic modulus 
20 distribution {E °} is obtained. 
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ly^bottom ^ 
Ur=«£/e-"^ ...(26) 

In the above Expressions, pi represents the external 
pressure vector on the body surface, and Sn represents the 
stress component orthogonal to the side face. The first 
5 Expression represents the boundary condition that the bottom 
is fixed, the second Expression represents the boundary 
condition that the stress distribution on the body surface 
matches the external pressure distribution, and the third 
Expression represents the boundary condition that each side 

10 face has a free end. Next, the displacement distribution 
{u"°} is substituted into the strain- displacement relation, 
whereby the strain distribution {e"°} corresponding to the 
elastic modulus distribution {E^} is obtained. Then, the 
strain distribution {e °} is substituted into the stress- 

15 strain relation, whereby the stress distribution {S °} 

corresponding to the elastic modulus distribution {E °} is 
obtained. 

Third, the elastic modulus distribution {E*^} is updated 
based upon the stress distribution calculated with the 
20 three-dimensional finite element method and the strain 
distribution {Ey} in the axial direction (y direction) 
estimated with the expanded combined autocorrelation method. 
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using the following Expression (27). 

...(21) 

Note that the above Expression is obtained by solving 
the aforementioned stress-strain relation for the Young's 
5 modulus E with 1=3=2. In the above Expression, k 
represents the iteration number. 

Fourth, the three-dimensional finite element analysis 
is consecutively made based upon the elastic modulus 
distribution thus updated and the aforementioned boundary 
10 conditions, whereby the stress distribution is updated. 

Then, the third and fourth processing are consecutively 
performed, whereby the elastic modulus distribution nears 
the actual elastic modulus distribution. Note that in the 
event that the elastic modulus distribution satisfies the 
15 following Expression (28), determination is made that 

convergence of the estimated elastic modulus distribution is 
achieved, and estimation processing ends. 




..(28) 



Here, 1 represents the element number, N represents the 
20 total number of the elements, and T represents the threshold 
value . 

Description has been made regarding the elastic modulus 
distribution reconstructing method using the three- 
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dimensional finite element model proposed in the present 
embodiment. With the present method, the elastic modulus 
distribution is estimated based upon the three-dimensional 
equilibrium equation. Accordingly, with the present method, 
5 calculation is made with assumptions which are closer to 
actual tissue than with the conventional methods, thereby 
enabling more precise estimation of the elastic modulus 
distribution. Furthermore, with the present method, the 
elastic modulus distribution is reconstructed based upon the 

10 strain distribution in the axial direction alone which can 
be estimated with high precision, thereby enabling 
reconstruction of the elastic modulus distribution in a sure 
manner. Note that with the present method, the three- 
dimensional distribution of the elastic modulus is estimated 

15 for the tissue, and accordingly, there is the need to employ 
a two-dimensional array ultrasonic probe, or there is the 
need to make mechanical scanning of a one -dimensional array 
ultrasonic probe in the slice direction, for making three- 
dimensional scanning of the tissue which is to be analyzed. 

20 Description will be made regarding the advantage of the 

reconstructing method for the elastic modulus distribution 
based upon the expanded combined autocorrelation method and 
the three-dimensional finite element model according to the 
present embodiment , which has been confirmed using 

25 simulation. Fig. 12 is a schematic diagram which shows the 
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procedure of the simulation. 

First, a tissue model having an elastic modulus 
distribution used for an estimation test Is created. Note 
that the tissue model contains scattering points which 
5 reflect ultrasonic echo signals therewlthln. Second, 

external pressure Is applied to the tissue model so as to 
cause compression thereof In computer simulation. Then, the 
displacement of each scattering point due to the compression 
Is calculated using the finite element method or the like. 

10 Third, the RF signals with and without application of 
pressure are simulated based upon the scattering-point 
displacement distribution with and without pressure applied 
to the tissue model. Fourth, the expanded combined 
autocorrelation method is applied to the simulated RF 

15 signals with and without application of pressure, whereby 
the tissue strain distribution is estimated. Fifth, the 
tissue elastic modulus distribution is estimated based upon 
the strain distribution estimated with the expanded combined 
autocorrelation method and the boundary conditions (external 

20 pressure distribution and so forth) determined for 

simulating compression of the tissue model using the elastic 
modulus distribution reconstructing method with the three- 
dimensional finite element model. 

While various elastic modulus distributions were used 

25 for the tissue models in this simulation, all cases of 
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elastic modulus distribution used here were assumed to 
exhibit isotropic elasticity. Note that the elastic modulus 
used in this simulation generally matches that of the 
mammary tissue corresponding to principal use of the tissue 
5 elastic modulus distribution measurement system according to 
the present embodiment. On the other hand, each tissue 
model contains scattering points therewithin for simulating 
the reflected RF signals with and without pressure applied 
to the tissue. Specifically, each tissue model contains the 

10 scattering points with average density of 500 points/cm^. 

Furthermore, the positions of the scattering points without 
application of pressure are determined using normal pseudo- 
random niimbers with an average of 1.0 and standard deviation 
of 0.3. Then, the positions of the scattering points under 

15 pressure are obtained by calculating the displacement of 

each scattering point without application of pressure based 
upon the results from the finite element analysis. Here, 
while the information with regard to the scattering points 
within the actual tissue is unknown, each parameter of the 

20 scattering point is determined such that the B-mode image 
created based upon the simulated RF signals is similar to 
the B-mode image of the actual tissue. 

With the present embodiment, the simulated RF signals 
with and without pressure applied to the tissue are 

25 calculated for each tissue model by convolution of the 
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scattering-point function with and without application of 
pressure and the point spread function of the ultrasonic 
system as represented by the following Expression (29). 

5 Here, ii(x, y, z) represents the RF signals without 

pressure applied to the tissue, iaCx, y, z) represents the 
RF signals under pressure applied to the tissue, h(x, y, z) 
represents the point spread function (impulse response) of 
the ultrasonic system, ti(x, y, z) represents the 

10 scattering-point function without pressure applied to the 
tissue, and t2(x, y, z) represents the scattering-point 
function under pressure applied to the tissue. Note that 
the scattering-point function represents the scattering 
coefficient thus predetermined at each scattering point , and 

15 represents a scattering coefficient of zero at positions 
other than the scattering points. On the other hand, the 
scattering-point function t2(x, y, z) under pressure applied 
to the tissue is calculated from the scattering-point 
function ti(x, y, z) without pressure applied to the tissue 

20 based upon the displacement of each scattering point due to 
strain of the tissue model. Note that the displacement of 
each scattering point due to compression of tissue is 
calculated by making linear interpolation of the 
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displacement vectors at the element nodes calculated with 
the finite element analysis. 

With the simulation according to the present embodiment, 
let us say that a non-focused ultrasonic system without 
5 damping of the ultrasonic wave is employed. That is to say, 
the point spread function h(x, y, z) is assumed to be 
constant for each point. Furthermore, let us say that the 
point spread function can be separated into functions for 
each direction as represented by the following Expression 
10 (30). 

Kx,y,z)^K{x)h^{y)Kiz) ^^^^ 

Here, hy(y) represents the point spread function in the 
direction of the ultrasonic beam. On the other hand, each 
of hx(x) and hz(z) represent the point spread function in the 

15 direction orthogonal to the ultrasonic-beam direction. 

Specifically, let us say that the x direction represents the 
in-plane direction (horizontal direction) of an ultrasonic 
tomographic image, and z direction represents the direction 
(slice direction) perpendicular to the ultrasonic 

20 tomographic image. Note that the point spread function in 
each direction is created based upon the reflected echo 
distribution obtained from actual measurement of the echo 
signals reflected from a wire target (a wire with a diameter 
of 0.13 mm extending in water) using an ultrasonic apparatus. 
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Fig. 13 shows diagrams for describing each point spread 
function used for simulation with an ultrasonic center 
frequency of 5.0 MHz. Fig. 13(A) shows a point spread 
function hy(y) in the axial direction. The point spread 
5 function hy(y) is obtained by multiplying the Gaussian 
function by a sine wave, and serves as an approximate 
distribution of the actual reflected echo distribution 
reflected from the wire target. On the other hand. Fig. 
13(B) shows the point spread function hx(x) in the 

10 horizontal direction, and Fig. 13(C) shows the point spread 
function hz(z) in the slice direction. Note that each of 
the aforementioned point spread functions are created using 
the Gaussian function so as to serve as an approximate 
distribution of the actual reflected echo distribution 

15 reflected from the wire target in the same way. The 

parameters of each function are varied corresponding to the 
center frequency, which will be described later in 
description of each simulation. 

Next, description will be made regarding the advantage 

20 of the expanded combined autocorrelation method as a 
displacement (strain) distribution estimating method 
according to the present embodiment, which has confirmed by 
simulation. First, description will be made regarding the 
advantage of estimating the displacement of the tissue in 

25 the horizontal direction, which is the advantage to the 
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combined autocorrelation method. 

Fig. 14 Is a schematic diagram which shows a tissue 
model. The tissue model Is formed with an outer size of 60 
mm X 60 mm (two-dimensional size), and with a uniform 
5 elastic modulus distribution. Then, compression of the 

tissue Is simulated so as to cause uniform strain of 3% In 
the axial direction. Let us say that a simple one- 
dlmenslonal elastic model Is employed as the tissue model In 
this simulation for evaluating only the advantage of the 

10 expanded combined autocorrelation method. Furthermore, 
compression of the tissue In the axial direction Is 
simulated with displacement of 0.0 mm to 1.4 mm In the 
horizontal direction for evaluating the advantage of 
handling displacement In the horizontal direction (relative 

15 displacement of the tissue In the horizontal direction as to 
the ultrasonic probe). Furthermore, let us say that simple 
parallel displacement of the tissue In the horizontal 
direction Is simulated, which corresponds to a situation 
wherein the ultrasonic probe has slipped as to the tissue. 

20 The RF signals are then simulated for the tissue model 

with and without strain of the tissue. The parameters used 
for simulation of the ultrasonic system Include: a center 
frequency of 5.0 MHz, a pulse width of 0.5 mm; a ultrasonic 
beam width of 1.0 mm; a scanning line pitch of 0.5 mm; and a 

25 sampling frequency of 30 MHz. Then, the strain distribution 



is estimated based upon the simulated RF signals with and 
without application of pressure with the expanded combined 
autocorrelation method proposed in the present embodiment. 
In addition, the strain distribution estimated with the 
5 combined autocorrelation method and the strain distribution 
estimated with the spatial correlation method were prepared 
for comparison. Note that a correlation window with the 
same size was applied to each method with the same search 
range, thereby allowing the user to make simple comparison 

10 of precision and so forth therebetween. Specifically, the 
expanded combined autocorrelation method and the spatial 
correlation method employ a two-dimensional window with a 
size of 1.6 mm (axial direction) x 2.5 mm (horizontal 
direction) for being applied in a two-dimensional search 

15 range of 5.6 mm (axial direction) x 7 . 5 mm (horizontal 
direction) . On the other hand, with the combined 
autocorrelation method for analyzing one -dimensional 
displacement, a one -dimensional correlation window with the 
same length of 1.6 mm in the axial direction was applied in 

20 the same range of 5.6 mm in the axial direction. 

Fig. 15 through Fig. 17 show estimation results of the 
strain distribution with each method.. Fig. 15 shows the 
error of the strain in the horizontal direction estimated 
with each method. Here, the reference symbol "O" 

25 represents the error with the combined autocorrelation 
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method, "□" represents the error with the expanded combined 
autocorrelation method, and "A" represents the error with 
the spatial correlation method. Note that the error of the 
estimated strain e is represented by the following 
5 Expression (31). 



Here, e ± represents the estimated strain, £i represents 
the actual strain (ideal value), i represents the element 
number, and N represents the total number of the elements. 

10 On the other hand. Fig. 16 shows the estimated strain 

distributions of the tissue containing displacement in the 
horizontal direction of 0.0 mm, using each method (combined 
autocorrelation method, expanded combined autocorrelation 
method, and spatial correlation method) . Fig. 17 shows the 

15 estimated strain distributions of the tissue containing 

displacement in the horizontal direction of 0.4 mm, using 
each method (combined autocorrelation method, expanded 
combined autocorrelation method, and spatial correlation 
method). Note that Fig. 16 and Fig. 17 show the average and 

20 the standard deviation of the estimated strain for each 
depth along the axial direction. 

As can be understood from the simulation results, with 
the conventional combined autocorrelation method (CA method). 



to 





i 



..(31) 
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relative displacement of the tissue in the horizontal 
direction exceeding the ultrasonic beam size (in this case, 
half the beam width, i.e., 0.5 mm) leads to rapid 
deterioration in estimation quality for the strain. On the 
5 other hand, it can be understood that the expanded combined 
autocorrelation method enables stable estimation of the 
strain regardless of the magnitude of the displacement in 
the horizontal direction. Furthermore, it can be understood 
that while the spatial correlation method enables stable 

10 estimation of the strain regardless of the magnitude of the 
displacement in the horizontal direction, as well, 
estimation results exhibit poor precision, i.e., exhibit 
twice or more the error as compared with the expanded 
combined autocorrelation method. Furthermore, making 

15 comparison between processing time of the aforementioned 

methods, while the expanded combined autocorrelation method 
requires processing time 3.6 times as great as with the 
combined autocorrelation method, the expanded combined 
autocorrelation method requires processing time only 1/(7.7) 

20 times as great as with the spatial correlation method, as 

shown in the following Table. As can be understood from the 
aforementioned results , the expanded combined 
autocorrelation method enables real-time calculation to a 
certain degree. 
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Method 


Processing time 


iMOinna.xizea 
processing time 


Combined 
autocorrelation 
method 


26 seconds 


1/(3.6) 


Expanded combined 
autocorrelation 
method 


1 minute 34 seconds 


1.0 


Spatial 
correlation 
method 


12 minutes 5 seconds 


7.7 



Next, description will be made regarding estimation 
results with application of pressure in the slant direction. 
Note that with the present estimation, simulation is made 
using a tissue model having a three-dimensional structure as 
with the actual tissue, unlike the aforementioned simulation 
using a simple tissue model having a two-dimensional uniform 
structure. Note that pressure should be applied to the 
tissue with an ultrasonic probe in the ultrasonic-beam 
direction (axial direction) in ideal measurement. Now, 
description will be made regarding the estimation results 
affected by compression of the tissue in the slant direction. 

Fig. 18 is a schematic diagram which shows a tissue 
model used for estimating the influence of compression of 
the tissue in a slant direction. As shown in Fig. 18(A), 
the tissue model has a three-dimensional structure with an 
outer size of 60 mm x 60 mm x 60 mm, and contains a 
cylindrical inclusion with a diameter of 15 mm and a length 
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of 60 mm, and with high rigidity. Let us say that the 
material surrounding the inclusion has an elastic modulus 
(Young's modulus) of 10 Kpa, and the inclusion has the 
elastic modulus three times as great as that of the 
5 aforementioned material, i.e., have an elastic modulus 

(Young's modulus) of 30 Kpa. Note that the aforementioned 
elastic moduli are determined based upon the elastic moduli 
of the mammary tissue and the elastic modulus of the mammary 
cancer, of which diagnosis is the principal object of the 

lb present invention. Then, compression of the tissue model is 
simulated in the following two situations. The one 
situation corresponds to compression of the tissue model by 
2% in the axial direction due to uniform external pressure 
of 200 Pa applied to the tissue model from the upper face 

15 along the axial direction as shown in Fig. 18(B). The other 
situation corresponds to compression of the tissue model in 
the slant direction due to uniform external pressure (200 Pa 
in the axial direction and 30 Pa in the horizontal 
direction) applied to the tissue model from the upper face 

20 along a slant direction as shown in Fig. 18(C). 

Then, the RF signals with and without application of 
pressure are simulated for the aforementioned two situations, 
following which the strain distribution is estimated with 
the expanded combined autocorrelation method. Note that the 

25 strain distribution is estimated with the combined 
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autocorrelation method and the spatial correlation method, 
as well, for comparison. Here, the correlation window 
having the same size Is applied to calculation In the same 
search range for each method, thereby allowing simple 
5 comparison. Note that the correlation window has the same 
size as with the simulation described above. Furthermore, 
the other parameters for simulating the RF signals are 
determined to be the same values as with the above -described 
simulation. That Is to say, the other parameters Include: a 
10 center frequency of 5.0 MHz, a pulse width of 0.5 mm; a 

ultrasonic beam width In the horizontal direction of 1.0 mm; 
a ultrasonic beam width in the slice direction of 2.0 mm; a 
scanning line pitch of 0.5 mm; and a sampling frequency of 
30 MHz. 

15 Fig. 19 and Fig. 20 show the simulation results In the 

aforementioned situations. Fig. 19 shows estimation results 
of the strain distribution In a simple situation wherein the 
tissue model Is compressed In the axial direction. Fig. 20 
shows estimation results of the strain distribution in a 

20 situation wherein the tissue model is compressed In a slant 
direction. Note that the strain distribution in the axial 
direction obtained by three-dimensional finite element 
analysis Is employed as the strain distribution estimated 
with ideal method for each situation, which serves as the 

25 actual strain distribution. Note that Fig. 19 and Fig. 20 
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are cross -sectional views of the tissue model taken along 
the center line thereof, which show estimation results- In 
Fig. 20, the strain distribution estimated with an ideal 
method exhibits left-right asymmetry. The reason is that 
5 pressure is applied in a slant direction. Specifically, in 
this case, the pressure is applied in a lower-right 
direction in the drawing, leading to large displacement in 
the horizontal direction at the upper-right part in the 
drawing . 

10 First, it has been confirmed that the expanded combined 

autocorrelation method exhibits generally the same 
estimation results for the strain distribution containing 
compression in the axial direction as with the combined 
autocorrelation method. On the other hand, in a situation 

15 wherein pressure is applied to the tissue in a slant 
direction, it has been also confirmed that while some 
regions cannot be estimated with the combined 
autocorrelation method due to great displacement in the 
horizontal direction, the expanded combined autocorrelation 

20 method enables stable estimation of the strain distribution 
regardless of the magnitude of the displacement in the 
horizontal direction as in the above description of the 
aforementioned simulation. On the other hand, it has been 
also confirmed that while the spatial correlation method 

25 enables stable estimation regardless of the magnitude of the 
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displacement in the horizontal direction, the spatial 
correlation method has significantly poorer estimation 
precision as compared with the expanded autocorrelation 
method. In addition to the simulation results described 
5 above, it has been confirmed that the expanded combined 
autocorrelation method has the advantage of estimation of 
strain distribution. 

As described above, it has been confirmed that the 
tissue- strain distribution can be estimated with high 

10 precision at high speed using the aforementioned expanded 
combined autocorrelation method. Now, description will be 
made regarding the estimation results obtained by simulation 
for the elastic modulus distribution reconstructing method 
with the three-dimensional finite element model proposed in 

15 the present embodiment. Note that the elastic modulus 
distribution reconstructing method is a method for 
estimating the elastic modulus distribution based upon the 
strain distribution, which serves as a method performed in 
the second stage in the tissue elasticity imaging system. 

20 The elastic modulus distribution reconstructing method 

proposed in the present embodiment has the principal 
function for estimating the elastic modulus distribution 
based upon the three-dimensional kinetic equilibrium 
equation. Now, the advantages of the elastic modulus 

25 distribution reconstructing method proposed in the present 
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embodiment are confirmed by making comparison between the 
elastic modulus distribution reconstructing method according 
to the present embodiment and the two-dimensional 
reconstructing method having the same processing except for 
5 calculation using the two-dimensional kinetic equilibrium 
equation. Note that calculation with the two-dimensional 
reconstructing method is made on the assumption that the 
plane strain occurs in the tissue. In the present 
simulation, two models each of which have a three - 

10 dimensional structure as with the actual tissue are employed 
as the tissue models as shown in Fig. 21. That is to say. 
Fig. 21 shows two tissue model examples each of which have a 
three-dimensional structure. Fig. 21(a) shows an inclusion- 
containing model containing an inclusion serving as a 

15 mammary tumor model. Specifically, the inclusion- containing 
model has an outer size of 100 mm x 100 mm x 100 mm, and 
contains a rigid inclusion with a diameter of 20 mm. Let us 
say that the inclusion has an elastic modulus of 30 kPa, and 
the material surrounding the inclusion has an elastic 

20 modulus of 10 kPa. The aforementioned elastic moduli are 
determined based upon the elastic modulus of the actual 
meumnary tissue in the same way as with the simulation 
described above. On the other hand, both the inclusion and 
the material surrounding the inclusion exhibit near- 

25 incompressibility , and accordingly, both are assumed to have 
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the same Polsson's ratio of 0.49. Fig. 21(b) shows a layer- 
structure model for simulating a layer- structure tissue such 
as a muscle. The layer- structure model has an outer size of 
100 mm X 100 mm x 100 mm, and contains a rigid layer with a 
5 thickness of 20 mm. Let us say that the rigid layer has an 
elastic modulus of 30 kPa, and the material surrounding the 
rigid layer has an elastic modulus of 10 kPa. Note that the 
layer- structure model has a uniform Poisson's ratio of 0.49, 
as well. 

10 Then, in a case of the inclusion-containing model shown 

in Fig. 21(a), compression under uniform external pressure 
of 100 Pa from the upper face of the model along the axial 
direction is simulated on a computer. On the other hand, in 
a case of the layer- structure model shown in Fig. 21(b), 

15 compression under uniform external pressure of 150 Pa from 
the upper face of the model along the axial direction is 
simulated on a computer 1 Note that compression under 
different external pressure is simulated for the 
aforementioned two models such that the saune strain of 

20 approximately 1% is simulated for each model. Then, the RF 
signals with and without application of pressure are 
simulated for each tissue model, and the strain distribution 
in the axial direction is estimated with the expanded 
combined autocorrelation method. Subsequently, the elastic 

25 modulus distribution is estimated based upon the estimated 



strain distribution in the axial direction and the boundary 
conditions determined for the simulation of compression, 
using the three-dimensional elastic modulus distribution 
reconstructing method. Also, the elastic modulus 
5 distribution is estimated based upon the same strain 

distribution in the axial direction and the same boundary 
conditions using the two-dimensional reconstructing method 
for comparison. Here, the parameters used for simulating 
the RF signals include: a center frequency of 3.7 5 MHz, a 

10 pulse width of 0.75 mm; a ultrasonic beam width in the 

horizontal direction of 2.0 mm; a ultrasonic beam width in 
the slice direction of 2.0 mm; and a scanning line pitch of 
2.0 mm. On the other hand, the parameters used for 
calculation with the expanded combined autocorrelation 

15 method include the size of the correlation window of 3.2 mm 
(in the axial direction) x 4 . 0 mm (in the horizontal 
direction), and the search range of 11.2 mm (in the axial 
direction) x 14.0 mm (in the horizontal direction). On the 
other hand, with the elastic modulus distribution 

20 reconstructing method using the three-dimensional finite 
element model, each tissue model is formed of 50,000 
rectangular- parallelepiped elements each of which have a 
size of 2 mm (in the axial direction) x 2 mm (in the 
horizontal direction) x 5 mm (in the slice direction). 

25 Fig. 22 through Fig. 25 show the simulation results. 
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Fig. 22 and Fig. 23 show the estimation results for the 
inclusion- containing model. On the other hand. Fig. 24 and 
Fig. 25 show the estimation results for the layer- structure 
model. Note that while the three-dimensional distribution 
5 of the elastic modulus can be estimated with the three- 
dimensional reconstructing method, each of Fig. 24 and Fig. 
25 show a cross-sectional view thereof taken along the 
center line of the model. The reason is that only the two- 
dimensional elastic modulus distribution can be estimated 

10 with the two-dimensional reconstructing method. Accordingly, 
Fig. 24 and Fig. 25 show a cross -sectional view of the 
estimation results taken along the center line of the model, 
thereby allowing the user to make comparison therebetween. 
On the other hand, estimation values for the three - 

15 dimensional reconstructing method and the two-dimensional 
reconstructing method obtained using each tissue model are 
listed in the following Table. 
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Error of 
elastic 
modulus in 

region 
surrounding 
core of 
model [ % ] 


Standard 
deviation 
in region 
surrounding 
core of 
model [ % ] 


Error of 
contrast 
in core 
of model 
[%] 


Inclusion- 
containing 
model 


Three- 
dimensional 
reconstructing 
method 


3.5 


15.5 


11.0 


Two- 
dimensional 
reconstructing 
method 


30.9 


17.9 


35.9 


Layer - 
structure 
model 


Three- 
dimensional 
r econs t rue t ing 
method 


8.5 


26.8 


3.1 


Two- 
dimensional 
reconstructing 
method 


24.9 


22.1 


43.5 



The estimation parameters used here include: the error of 
elastic modulus es in the region surrounding the core of the 
model; the standard deviation SDs in the region surrounding 
the core of the model; and the error of the contrast ec in 
the core region of the model, which are defined by the 
following Expressions . 



SJD, 



.(32) 
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Note that in the aforementioned Expressions, S 
represents the sum in the region surrounding the core, E 
represents the estimated elastic modulus, E represents the 
actual elastic modulus, Ns represents the total number of 
5 the elements in the region surrounding the core, E"s 

represents the average of the elastic modulus in the region 
surrounding the core of the model, E c represents the 
estimated elastic modulus in the core region of the model, 
Ec represents the actual elastic modulus in the core region 

10 of the model, and Es represents the actual elastic modulus 
in the region surrounding the core of the model. 

As can be understood from the aforementioned simulation 
results, it has been confirmed that the elastic modulus 
distribution reconstructing method with the three- 

15 dimensional finite element model proposed in the present 
embodiment has the advantage of estimating the elastic 
modulus distribution with higher precision than with the 
two-dimensional reconstructing method formed on the 
assumption that plane stress occurs within the tissue. 

20 While three-dimensional reconstructing method enables 

precise estimation of the elastic modulus distribution, the 
elastic modulus distribution estimated with the two- 
dimensional reconstructing method exhibits smaller values 
than those of the actual elastic modulus distribution. It 

25 is needless to say that the reason is that with the two- 
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dimensional reconstructing method, calculation In the 
direction orthogonal to the two-dimensional calculation 
plane Is limited. In the present evaluation. It has been 
clearly confirmed that the elastic modulus distribution 
5 reconstructing method with the three-dimensional calculation 
corresponding to the actual tissue Is suitably employed for 
analyzing actual tissue. 

Next, description will be made regarding a specific 
configuration of an ultrasonic diagnosis system employing 

10 the aforementioned expanded combined autocorrelation method 
and the aforementioned elastic modulus reconstructing method 
with the three-dimensional finite element model. Fig. 26 Is 
a diagram which shows a basic configuration of the 
ultrasonic diagnosis system. The ultrasonic diagnosis 

15 system comprises a three-dimensional ultrasonic scanner 281, 
an ultrasonic diagnosis apparatus 282, a personal computer 
283, a pulse motor controller 284, a pulse motor 285, a 
pressure gauge 286, and the like. The three-dimensional 
ultrasonic scanner 281 has functions for casting ultrasonic 

20 pulse signals toward the tissue, and for receiving 

ultrasonic echo signals reflected from the tissue. Note 
that the system employs the elastic modulus reconstructing 
method with the three-dimensional finite element model, and 
accordingly, the system requires three-dimensional data 

25 within the tissue. Accordingly, the ultrasonic diagnosis 
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system has a configuration for making three-dimensional 
scanning with the three-dimensional ultrasonic scanner 281. 
The ultrasonic diagnosis apparatus 282 has functions for 
controlling the three-dimensional ultrasonic scanner 281, 
5 and for displaying real-time ultrasonic B-mode images, 

thereby allowing the user to determine the position of the 
region which is to be measured. Note that a conventional 
ultrasonic diagnosis apparatus may be employed as the 
ultrasonic diagnosis apparatus 282 without modification. 

10 The ultrasonic diagnosis apparatus includes a full-digital 
device having frame memory for temporarily storing the 
measured RF signals. The personal computer 283 receives the 
RF signals measured by the ultrasonic diagnosis apparatus 
282, performs processing (processing proposed in the present 

15 embodiment described above) for estimating the tissue 

elastic properties, and displays the estimation results. 
The pulse motor 285 has a function for controlling the 
pressure applied to the tissue. The pulse motor 285 is 
fixed at the tip of a stationary arm, and the three- 

20 dimensional ultrasonic scanner 281 is mounted on the moving 
portion of the pulse motor 285. The system has a 
configuration wherein the pulse motor controller 284 
controls the pulse motor 285 for adjusting the position of 
the ultrasonic scanner 281 in the vertical direction, 

25 thereby allowing the user to adjusting the pressure applied 
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to the tissue which causes small compression thereof by 
several percent with high precision. The pressure gauge 286 
has a function for measuring the pressure on the body 
surface; the pressure serving as the boundary condition 
5 required for reconstructing the elastic modulus 

reconstruction. Note that the pressure gauge 286 is 
positioned between the ultrasonic scanner 281 and the body 
surface. Here, the pressure measured with the pressure 
gauge 286 is used as the pressure on the body surface on the 

10 assumption that compression of the tissue through the 
ultrasonic scanner 281 causes the uniform pressure 
distribution on the body surface. 

Fig. 27 is a diagram which shows a specific 
configuration of the ultrasonic scanner 281 included in the 

15 ultrasonic diagnosis system. The three-dimensional 

ultrasonic scanner 281 does not have a two-dimensional array 
configuration wherein ultrasonic transducers are two- 
dimensionally arrayed, but has a three-dimensional scanning 
configuration wherein mechanical scanning of a two- 

20 dimensional convex scanning probe is made in the slice 
direction in water. 

The ultrasonic diagnosis system shown in Fig. 26 is 
designed principally for diagnosis of mammary cancer, and 
accordingly, the property parameters of the ultrasonic 

25 scanner are determined for diagnosis of mammary cancer. 



1 



Specifically, the principal property parameters of the 
ultrasonic scanner used here include: the ultrasonic center 
frequency of 7.5 MHz, the sampling frequency of 30 MHz, the 
number of scanning lines of 71, the nximber of frames of 44, 
5 the transducer scanning angle of 30^ , and the period of time 
for a three-dimensional scanning cycle of 0.5 seconds. Note 
that scanning of the convex probe is made in the slice 
direction by the aforementioned scanning angle of the 
transducer- Also, the aforementioned number of freunes means 

10 the number of the scanning images (frames) acquired for a 

single cycle of scanning of the convex probe. On the other 
hand, the properties of the ultrasonic pulse obtained by 
actual measurement using a wire target in water include: the 
pulse width of approximately 0.5 mm; the beam width in the 

15 horizontal direction of approximately 1.5 mm; and the beam 
width in the slice direction of approximately 2.6 mm. 

Next, description will be made regarding the operation 
of the ultrasonic diagnosis system shown in Fig. 26, which 
makes measurement of the elastic properties. First, the 

20 user moves the three-dimensional scanner 281 mounted on the 
arm so as to make positioning of the ultrasonic scanner 281 
at a desired portion which is to be measured while 
monitoring the real-time B-mode images obtained by the 
ultrasonic diagnosis apparatus 282. Note that at the time 

25 of positioning of the ultrasonic scanner 281, three- 
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dimensional scanning of the ultrasonic scanner 281 Is not 
made (I.e., mechanical scanning of the convex probe Is not 
made), and only the B-mode Image corresponding to the center 
plane of the scanner Is displayed on the ultrasonic 
5 diagnosis apparatus 282. The reason Is that with the 

ultrasonic diagnosis apparatus 282 used here, the real- time 
B-mode Images cannot be displayed at the time of three- 
dimensional scanning. It Is needless to say that an 
ultrasonic diagnosis apparatus may be employed, which has a 

10 function wherein the real-time B-mode Images can be 

displayed at the time of three-dimensional scanning, thereby 
allowing the user to make positioning of the three- 
dimensional ultrasonic scanner 281 while making three- 
dimensional scanning. Following positioning of the 

15 ultrasonic scanner 281 at a desired portion which Is to be 
measured, the user fixes the moving portion of the arm in 
order to fix the ultrasonic scanner 281. Then, the system 
makes measurement of the three-dimensional RF signals 
without pressure applied to the tissue. Note that upon the 

20 user pressing the button for three-dimensional scanning, 

three-dimensional scanning is automatically made. Here, the 
period of time for a single cycle of the three-dimensional 
scanning is only 0 . 5 seconds . The measured RF data without 
application of pressure is stored in the frame memory within 

25 the ultrasonic apparatus. Subsequently, upon the user 
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pressing the button one time for operating the pulse motor 
controller 284 for controlling compression of the tissue, 
the pulse motor 285 fixed to the arm moves the ultrasonic 
scanner 281 by a predetermined distance, thereby compressing 
5 the tissue through the ultrasonic scanner 281. Subsequently, 
the motor 285 stops while maintaining compression of the 
tissue. In this state, the user presses the button again 
for three-dimensional scanning, whereby the RF data under 
pressure applied to the tissue is measured. Note that the 

10 RF data under pressure applied to the tissue is stored in 

the frame memory included in the ultrasonic apparatus 282 in 
the same way as with the RF data without application of 
pressure. Also, the pressure applied to the tissue is 
measured with the pressure gauge mounted at the end of the 

15 ultrasonic scanner 281. Then, measurement ends, and the 

pressure applied to the tissue is released, following which 
the subject is freed. 

Following freeing of the subject, the system accesses 
the frame memory included in the ultrasonic diagnosis 

20 apparatus 282 through the personal computer 283, and the RF 
data with and without pressure applied to the tissue is 
stored on a hard disk included in the personal computer 283. 
Such processing is performed since the frame memory included 
in the ultrasonic apparatus has a function for temporarily 

25 storing the data, i.e., has only a capacity for storing the 
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data for a single measurement cycle. The personal computer 
283 executes a program for calculation using the expanded 
combined autocorrelation method and the elastic modulus 
distribution reconstructing method with the three- 
5 dimensional finite element model in order to estimate the 
strain distribution and the elastic modulus distribution 
based upon the RF data with and without pressure applied to 
the tissue. Then, the personal computer 283 displays the B- 
mode image, the strain image, and the elastic modulus image, 

10 on a monitor, by executing a display program. 

The ultrasonic diagnosis system using the strain 
distribution display method and the elastic modulus 
distribution display method according to the present 
invention has the advantage of enabling reconstruction of 

15 the elastic modulus distribution based upon the strain 
distribution in the ultrasonic beam direction (axial 
direction) alone, as well as the advantage of enabling 
estimation of the displacement distribution regardless of 
displacement in the horizontal direction. 

20 Note that while description has been made regarding an 

arrangement wherein the envelope signals are employed, the 
present invention is not restricted to the aforementioned 
arrangement, rather, any parameter which represents the 
relation between the wave properties including the amplitude, 

25 the wave height, and the number of waves, may be employed. 



